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This  material  was  originally  published  in  December  1988  as  a 
dissertation,  in  partial  fulfillment  of  the  requirements  for  the  degree  of  Doctor  of 
Philosophy  in  Electrical  Engineering,  The  University  of  Texas  at  Austin. 
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CHAPTER  1 


INTRODUCTION 


In  the  field  of  underwater  sound,  the  ability  to  model  acoustic  propagation 
is  an  essential  part  of  most  practical  applications.  For  example,  the  development 
of  passive  and  active  sonar  systems,  the  estimation  of  ranges  at  which  targets  can 
be  detected,  and  the  study  of  the  effects  of  ambient  noise  sources  all  require  an 
understanding  of  how  sound  propagates  in  complex  ocean  environments. 

In  recent  years,  the  advent  of  ever  more  powerful  computers  has  enabled 
researchers  to  develop  and  implement  sophisticated  propagation  models  that  nu¬ 
merically  solve  the  acoustic  wave  equation  or  approximations  to  it.  Such  models 
rely  on  a  thorough  description  of  the  environment,  which  must  include  the  geome¬ 
try  (bathymetry,  layering  in  the  ocean  bottom,  and  sea  surface  state,  for  example), 
as  well  as  the  acoustic  properties  of  the  media  (sound  speed,  density,  and  attenua¬ 
tion,  for  example).  However,  even  if  the  environment  were  known  exactly,  models 
typically  have  limitations  associated  with  them,  such  as  the  ability  to  properly 
account  for  (1)  dependence  of  the  environmental  parameters  on  range,  (2)  propa¬ 
gation  in  three-dimensional  environments,  (3)  energy  propagating  at  angles  near 
vertical,  (4)  shear  wave  propagation  in  the  ocean  bottom,  (5)  scattering  from 
rough  surfaces,  and  (6)  nonlinear  effects.  Practical  limitations,  such  as  the  range 
of  frequencies  that  can  be  modeled,  often  arise  simply  because  the  computation 
time  is  too  large  or  the  storage  capaci'y  of  the  computer  is  exceeded. 

A  problem  with  the  development  of  any  model  is  the  evaluation  of  its 
accuracy.  There  are  always  factors  that  can  introduce  error  in  the  results  of  a 
propagation  model,  such  as  approximations  in  the  theoretical  formulation,  errors 
arising  from  numerical  methods  used  in  the  computer  algorithm,  or  simply  er¬ 
rors  in  the  coding  of  the  program.  One  way  to  verify  the  accuracy  of  a  model 
is  to  compare  its  solutions  with  those  of  another  model  that  relies  on  a  different 
method.  This  technique  for  model  verification  requires  the  definition  of  certain 
“benchmark”  problems  to  which  different  models  can  be  applied  and  their  solu¬ 
tions  compared.  Recently,  the  Underwater  Acoustics  community  of  the  Acoustical 
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Society  of  America  proposed  benchmark  problems  that  address  propagation  in 
range-dependent  environments.  More  specifically,  the  problems  consist  of  a  point 
source  in  a  two-dimensional,  isovelocity  wedge  of  shallow  water,  surrounded  above 
by  air  and  below  by  sand.  One  concrete  result  of  the  present  research  is  a  new, 
accurate,  and  physically  intuitive  method  for  acoustic  propagation  modeling  that 
can  be  applied  to  these  benchmark  problems. 

The  goal  of  the  research  in  this  dissertation  has  been  to  develop  a  method 
based  on  ray  theory  that  is  accurate  in  shallow  water  ocean  environments.  Since 
classical  ray  theory  is  based  on  a  high-frequency  approximation  to  the  wave  equa¬ 
tion,  it  is  not  valid  in  shallow  water,  where  the  water  depth  is  typically  on  the 
order  of  tens  of  wavelengths  or  less.  In  order  to  adapt  classical  ray  theory,  we 
begin  with  an  exact  formulation  for  the  reflection  and  transmission  of  a  spheri¬ 
cal  wave  incident  on  a  single  plane  interface.  Preserving  as  much  as  possible  the 
intuitive  notions  of  ray  theory,  we  develop  a  method  for  solving  the  single  inter¬ 
face  problem  accurately,  even  it  low  frequencies.  We  then  extend  the  method 
to  calculate  the  field  due  to  rays  that  undergo  multiple  reflections.  Using  these 
tools,  we  can  calculate  the  total  field  in  ocean  environments  with  flat  or  sloping 
bottoms  by  summing  together  all  of  the  significant  rays.  Unlike  most  numerical 
models,  which  provide  little  or  no  physical  insight  into  the  paths  and  mechanisms 
of  propagation,  the  ray  method  developed  here  can  be  analyzed  in  terms  of  ray 
arrival  angles,  travel  times,  and  propagation  paths. 

This  dissertation  is  organized  as  follows.  In  Chapter  2  the  basic  method 
for  finding  the  reflected  and  transmitted  fields  due  to  a  point  source  in  the  presence 
of  a  single  plane  interface  is  developed.  The  extension  to  oceanic  environments 
with  flat  and  sloping  bottoms  is  made  in  Chapter  3.  In  Chapter  4  we  address  in 
detail  the  benchmark  wedge  problems  and  show  comparisons  between  solutions 
obtained  from  the  ray  model  and  from  a  two-way  coupled  mode  model.  The 
problem  of  propagation  in  a  three-dimensional  wedge,  where  sound  travels  across 
the  sloping  bottom  rather  than  straight  up  or  down  the  slope,  is  discussed  in 
Chapter  5.  In  Chapter  6,  we  demonstrate  how  the  methods  we  have  developed 
for  modeling  a  point  source  can  be  used  to  model  a  directional  source.  Several 
examples  of  the  reflected  and  transmitted  fields  due  to  a  beam  incident  on  a  fluid- 
fluid  interface  near  the  critical  angle  are  given.  Concluding  remarks  are  made  in 
Chapter  7. 

Three  of  the  chapters  of  this  dissertation  are  based  on  manuscripts  that 
were  submitted  for  publication  in  the  Journal  of  the  Acoustical  Society  of  Amer¬ 
ica  in  1988.  Chapter  2  is  based  on  the  manuscript'  entitled  “Complex  ray  meth¬ 
ods  for  acoustic  interaction  at  a  fluid  fluid  interface,”  Chapter  3  is  based  on  i  he 
manuscript2  entitled  “Hay  methods  for  flat  and  sloping  oceanic  waveguides,”  and 
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Chapter  4  is  based  on  the  manuscript3  entitled  “Ray  model  solutions  to  the  bench¬ 
mark  wedge  problems.”  The  results  in  the  last  manuscript  were  first  reported  at 
the  115th  Meeting  of  the  Acoustical  Society  of  America  in  Seattle,  Washington, 
in  May  1988. 4 


CHAPTER  2 


COMPLEX  RAY  METHODS  FOR  ACOUSTIC  INTERACTION  AT 

A  FLUID-FLUID  INTERFACE 


2-1  INTRODUCTION 

One  of  the  most  basic  problems  in  acoustics  is  that  of  the  reflection  and 
transmission  of  a  spherical  wave  at  a  plane,  fluid-fluid  interface.  In  the  case  of  un¬ 
derwater  acoustics,  this  problem  represents  an  idealized  model  for  the  interaction 
of  a  point-source  field  with  the  ocean  bottom.  More  specifically,  one  isovelocity, 
constant- density  halfspace  (water)  overlies  another  (sand,  for  example).  In  shal¬ 
low  water  ocean  environments,  the  lower  medium  usually  has  the  higher  velocity, 
which  results  in  the  presence  of  a  critical  angle  in  the  plane  wave  reflection  coef¬ 
ficient.  Plane  waves  incident  on  the  interface  at  angles  “steeper”  than  the  critical 
angle  (closer  to  normal  incidence)  are  partially  reflected  back  into  the  water  and 
partially  transmitted  into  the  bottom.  Plane  waves  incident  on  the  interface  at 
angles  “shallower”  than  the  critical  angle  (closer  to  grazing  incidence)  are  totally 
reflected  back  into  the  water.  When  the  incident  wave  is  a  spherical  wave,  en¬ 
ergy  incident  near  the  critical  angle  produces  a  wave  that  travels  parallel  to  the 
interface  in  the  bottom.  As  it  propagates,  it  continually  excites  a  wave  called  the 
lateral  wave  or  head  wave ,  which  travels  at  the  critical  angle  in  the  water. 

The  initial  step  in  the  mathematical  analysis  of  the  two-fluid  model  is 
to  represent  the  incident,  spherical  wave  as  a  sum  of  (integral  over)  plane  waves. 
An  integral  expression  for  the  reflected  or  transmitted  field  at  a  given  receiver 
point  is  then  obtained  by  multiplying  each  plane  wave  in  the  integrand  by  the 
appropriate  plane  wave  reflection  or  transmission  coefficient.  In  the  classic  work 
of  Brekhovskikh,5  the  author  primarily  uses  the  method  of  steepest  descent  to 
solve  the  integral.  Briefly  stated,  this  method  involves  finding  the  saddle  point, 
deforming  the  contour  of  integration  into  the  steepest  descent  path,  and  deriving 
an  approximate  solution  to  the  integral.  When  viewed  in  terms  of  rays,  the  saddle 
point  corresponds  to  the  angle  of  the  eigenray,  a  specific  ray  connecting  the  source 
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and  receiver.  The  major  complication  in  this  analysis  arises  from  the  presence  of 
the  critical  angle,  which  in  mathematical  terms  corresponds  to  a  branch  point 
singularity  in  the  reflection  or  transmission  coefficient.  When  the  saddle  point 
lies  sufficiently  close  to  the  critical  angle,  the  usual  methods  for  estimating  the 
field  break  down. 

Various  improvements  to  the  solution  in  the  neighborhood  of  the  critical 
angle  have  been  developed.  In  Sec.  31  of  Ref.  5,  Brekhovskikh  derives  expres¬ 
sions  for  the  reflected  and  lateral  wave  fields  in  the  vicinity  of  the  critical  angle. 
Bleistein6  derives  uniform  asymptotic  expansions  for  integrals  involving  station¬ 
ary  points  near  algebraic  singularities.  In  the  present  case,  a  saddle  point  lies 
near  a  branch  point  of  order  1/2.  Stickler7  refines  the  formulation  by  allowing 
for  the  possibility  that  the  critical  angle  is  close  to  grazing  incidence.  Chin-Bing 
and  Davis8  apply  Bleistein’s  technique  to  the  lateral  wave  expression  and  point 
out  errors  in  Brekhovskikh’s  treatment.  In  all  of  these  investigations,  asymptotic 
expansions  involving  parabolic  cylinder  functions  are  derived. 

In  the  present  work,  we  wish  to  develop  a  method  for  finding  the  reflected 
and  transmitted  fields  that  (1)  can  be  interpreted  in  terms  of  rays,  (2)  is  applicable 
even  at  very  low  frequencies,  and  (3)  can  be  extended  to  environments  where  more 
than  one  interface  is  present.  The  above-mentioned  asymptotic  expansions  do  not 
satisfy  these  criteria  because  (1)  there  is  no  physical  interpretation  of  the  terms 
in  the  expansion,  (2)  an  expansion  with  a  fixed  number  of  terms  will  fail  at  a 
sufficiently  low  frequency,  and  (3)  the  treatment  of  multiply  reflected  fields  would 
require  a  rederivation  of  the  asymptotic  expansion. 

The  approach  we  take  in  this  dissertation  is  based  on  the  saddle  point 
method  outlined  by  Brekhovskikh.5  In  one  variation  on  this  method,  Brekhovskikh 
(Sec.  31.3)  allows  the  phase  of  the  reflection  coefficient  to  influence  the  location 
of  the  saddle  points.  The  result  is  a  phenomenon  commonly  referred  to  as  beam 
displacement,  where  totally  reflected  rays  have  a  horizontal  displacement  along 
the  interface  before  they  are  reflected  back  into  the  water.  (Since  the  term  “beam 
displacement”  originated  from  the  displacement  of  a  beam  incident  on  an  interface, 
a  more  appropriate  term  for  the  present  case  would  be  ray  displacement.)  In  recent 
years,  the  phenomenon  of  beam  displacement  has  been  included  in  ray  theory 
calculations,  and  the  resulting  method  successfully  applied  to  the  flat,  isovelocity, 
penetrable-bottom  waveguide  (the  Pekeris  model)9,10,11  and  to  the  sloping-bottom 
waveguide.12 

The  one  deficiency  in  the  beam  displacement  formulat  ion,  however,  is  that 
it  is  not  accurate  in  a  certain  “transition”  region,  where  eigenray  angles  are  close 
to  the  critical  angle.  At  shorter  ranges,  the  field  can  be  expressed  in  terms  of  a 
single  partially  reflected  eigenray,  and  at  longer  ranges,  the  field  can  be  expressed 
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as  a  sum  of  two  separate  rays:  a  totally  reflected  eigenray  and  a  lateral  wave 
eigenray.  But  in  the  transition  region,  the  field  is  more  complicated  and  cannot  be 
expressed  a s  a  simple  ray  field  or  even  as  a  sum  of  simple  ray  fields.  In  this  region 
the  ray  formulation  breaks  down  —  it  predicts  infinite  intensity  at  caustics,  points 
in  space  where  two  rays  coalesce  into  one.  Special  caustic  correction  factors11,12 
based  on  uniform  asymptotics  are  used  to  attempt  to  “patch  up”  the  ray  solution 
in  this  problematic  region. 

In  the  approach  developed  in  this  chapter,  the  key  modification  is  to  al¬ 
low  the  entire  reflection  coefficient  (amplitude  as  well  as  phase)  to  influence  the 
locations  of  the  saddle  points.  The  result  is  that  the  saddle  points  become,  in  gen¬ 
eral,  complex,  corresponding  to  eigenrays  that  have  complex  angles  of  incidence  0. 
In  the  neighborhood  of  the  critical  angle,  the  saddle  point  approximation  to  the 
integral  still  fails.  However,  in  this  region  we  find  the  exact  field  by  numerically 
integrating  along  the  steepest  descent  path  in  the  complex  0-plane.  The  num¬ 
ber  of  integration  points  is  made  independent  of  frequency  by  choosing  step  sizes 
appropriately.  In  work  done  concurrently  with  the  present  research,  Plumpton 
and  Tindle13  formulate  an  identical  complex  saddle  point  criterion  for  the  single 
reflection  problem  and  develop  improved  caustic  correction  procedures. 

The  notion  of  incorporating  the  reflection  and  transmission  coefficients 
into  the  saddle  point  formulation  is  attractive  because  it  allows  the  full  effects 
of  the  interaction  with  the  interface  to  be  included.  For  example,  in  the  case 
of  the  water-sand  interface,  our  method  accounts  for  the  reflected  lateral  wave 
field  as  well  as  the  evanescent  wave  field  in  the  bottom.  In  cases  where  the  total 
field  is  composed  of  several  different  ray  fields,  the  individual  contributions  can 
be  isolated. 

In  this  chapter  we  develop  our  ray  method  while  solving  for  the  reflected 
field,  and  then  apply  it  in  a  straightforward  manner  to  the  transmitted  field.  In 
Sec.  2-2  we  derive  the  plane  wave  integral  for  the  reflected  field.  In  Sec.  2-3  we 
review  the  method  of  steepest  descent  and  apply  it  in  the  usual  way  to  the  reflected 
field  integral  to  obtain  classical  ray  theory.  In  Sec.  2-4  we  obtain  the  complex 
ray  formulation  by  including  the  plane  wave  reflection  coefficient  in  the  saddle 
point  criterion.  The  resulting  saddle  point  structure  is  discussed  in  Sec.  2-5,  and 
the  technique  for  evaluating  the  field  is  summarized  in  Sec.  2-6.  In  Sec.  2-7  we 
interpret  the  saddle  point  formulation  in  terms  of  eigenrays.  Section  2-8  contains 
the  solution  for  and  analysis  of  the  transmitted  field.  Concluding  remarks  are 
made  in  Sec.  2-9. 
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2-2  PLANE  WAVE  INTEGRAL  REPRESENTATION  OF  THE  RE¬ 
FLECTED  FIELD 


In  this  section  we  derive  the  plane  wave  integral  representation  for  the 
reflected  field.  The  derivation  follows  the  treatment  by  Brekhovskikh  in  Secs.  26 
and  28  of  Ref.  5.  The  first  step  is  to  express  the  pressure  field  from  an  acoustic 
point  source  in  free  space  as  an  integral  over  plane  waves.  This  type  of  formu¬ 
lation  is  attractive  because  the  interaction  of  individual  plane  waves  with  planar 
interfaces  is  mathematically  tractable  and  well  understood. 

We  begin  with  the  linear,  scalar  wave  equation  for  the  acoustic  pressure 
p  in  a  homogeneous  medium  of  sound  speed  c: 


1  d2p 

72  dF 


(2-1) 


We  consider  the  following  solution  to  Eq.  (2-1),  which  represents  the  field  due  to 
a  harmonic  point  source  at  the  origin  of  a  Cartesian  coordinate  system: 


i(kR-uit)  - - 

p{x,y,z)  =  p0 - - -  ,  where  R  =  yjx2  +  y2  +  z2  f  0  .  (2-2) 


The  quantity  p0  in  Eq.  (2-2)  is  a  constant  (in  units  pressure-distance)  that  gives 
the  pressure  at  unit  distance  from  the  origin.  After  substituting  Eq.  (2-2)  into 
Eq.  (2-1),  carrying  out  the  partial  derivatives  with  respect  to  t,  and  expressing 
the  spatial  derivatives  in  Cartesian  coordinates,  we  obtain  the  scalar  Helmholtz 
equation 


dx 2  dy2 


k2  p  =  0  , 


(2-3) 


where  the  wavenumber  is  defined  as  k  =  ujc. 

In  order  to  derive  a  plane  wave  representation  of  the  spherical  wave  in 
Eq.  (2-2),  we  first  write  the  pressure  field  in  the  z  =  0  plane  as 


ikr 

p(*,y)  =  — 

r 


where  r  =  yjx2  +  y2  ^  0 


(2-t) 


In  Eq.  (2-4),  the  e-,wt  time  dependence  of  the  field  has  been  suppressed,  and,  for 
conveniv  ice,  we  have  taken  p0  to  be  unity.  We  now  express  the  field  in  the  ( x,y ) 
plane  as  a  double  Fourier  integral  over  the  transform  variables  ( kx,kv ): 


e 


ikr 


r 


oo 

J j  A(kr,ky)e'lk'x+k»v'dkxdky 

—  OO 


(r?  0) 


(2-5) 
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where 


*<*.,*,)  ={^lJ  v  di d« 


(2-6) 


The  integral  in  Eq.  (2-6)  is  solved  by  converting  from  (x,  y)  coordinates 
to  polar  coordinates  (r,  <f>)  and  solving  the  two  integrals  analytically  (see  Ref.  5, 
pp.  228-229).  The  result  is 


A(kx,ky)  = 


2 xJk*  -k*-k* 


(2  7) 


Using  Eqs.  (2-5)  and  (2-7)  in  Eq.  (2-4),  we  may  write  the  pressure  field  in  the 
2  =  0  plane  as 


OC 

p(*.y)  =  // 


iL  V*a  -  *2  -  kl 


c>(k,r+kvy)  dki  dky  (r^O)  .  (2  8) 


We  now  construct  an  expression  for  the  spherical  wave  field  off  the  2  =  0 
plane  by  formally  multiplying  the  integrand  of  Eq.  (2-8)  by  e±,k‘*  for  ±z  >  0: 


p(x,  y,  z)  =  /  /  - dk  dk  z  >  o 

JJ  27r./it2  -  kl  —  kl  y 

-oo  V  1  V 


(2-9) 


This  “continuation”  is  valid  because  the  integrand  of  Eq.  (2-9),  and  therefore  p 
itself,  satisfies  the  Helmholtz  equation  (2-3)  and  reduces  to  Eq.  (2-8)  when  2  =  0. 
Substitution  of  Eq.  (2-9)  into  the  Helmholtz  equation  yields  the  relation 


kl  +  k2y  +  kl  =  k2  . 

Recognizing  the  radical  in  Eq.  (2-  9)  to  be 

=  vA2  -  kl  -kl  , 

we  may  rewrite  the  expression  for  the  spherical  wave 


(2-10) 


(2-11) 


exkR  f  7  i 

P(x,j/,2)  =  —  =  jj  —  c*k'r+k”'+k'*'dktdky  2>0  .  (2  12) 


Die  exponential  term  in  Eq.  (2  12)  represents  a  plain*  wave  whose  di¬ 
rection  is  given  by  the  three  wave  numbers  kr,  ky ,  and  kr,  the  last  of  which  is 
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dependent  on  the  first  two  by  way  of  Eq.  (2-11)  (see  Fig.  2-1).  Each  point  in 
the  (kx,ky)  plane  can  be  viewed  as  representing  a  plane  wave  that  contributes  to 
the  spherical  wave  representation  in  Eq.  (2-12).  The  fraction  in  the  integrand 
of  Eq.  (2-12)  can  be  viewed  as  a  weighting  function  for  each  plane  wave.  Points 
inside  the  circle  of  radius  k  in  the  (kx,ky)  plane  correspond  to  the  usual,  propa¬ 
gating  (homogeneous)  plane  waves  since,  by  Eq.  (2-11),  kz  is  real.  Points  outside 
the  circle  correspond  to  “inhomogeneous”  plane  waves,  for  which  k2  is  imagi¬ 
nary.  By  substituting  a  positive  imaginary  value  for  kz  into  the  plane  wave  form 
et(k,x+kyy+kzz) ,  can  ke  geen  ^at  inhomogeneous  plane  waves  are  attenuated  in 
the  positive  2-direction.  The  inclusion  of  the  inhomogeneous  plane  waves  is  nec¬ 
essary  to  mathematically  complete  the  description  of  the  spherical  wave  in  terms 
of  plane  waves. 


Figure  2-1  Change  of  variable  from  rectangular  coordinates  (kI<ky,k[)  to  polar 

COORDINATES  (k,0,<t>).  THE  DIFFERENTIAL  ELEMENT  OF  AREA  IN  THE  (kx,ky)  PLANE  IS 
SHOWN. 


We  now  convert  the  wave  number  coordinates  from  rectangular  ( kx ,  ky,  k2) 
to  polar  (k,0,<f>)  by  way  of  the  relations 

kx  =  k  sin  6  cos  (f> 

ky  =  k  sin  0  sin  <f>  .  (2-13) 

kz  =  k  cos  0 
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These  relations  are  shown  graphically  in  Fig.  2-1.  From  Fig.  2-1,  the  differential 
element  of  area.  dkx  dkv  transforms  to 

dkx  dky  «-»  (&sin  0 d4>)  (kcos  Odd)  ,  (2  14) 

where  the  first  term  reflects  the  ^-dependence  and  the  second  term  the  ^-depend¬ 
ence.  To  cover  the  whole  ( kx ,  ky)  plane,  <f>  must  vary  from  0  to  2tt.  The  inside  of 
the  circle  of  radius  k  is  covered  when  0  varies  from  0  to  7r/2,  while  the  outside  of 
the  circle  is  covered  when  0  varies  from  tt/2  to  7r/2  —  zoo.  The  imaginary  part  of  0  is 
chosen  to  go  from  0  to  — t oo  in  order  to  make  the  amplitudes  of  the  inhomogeneous 
plane  waves  in  the  integrand  of  Eq.  (2-12)  exponentially  decreasing  rather  than 
increasing,  which  can  be  seen  by  writing 


k;  =  k  cos  6  =  k  cos  (On  +  i0\)  =  A;(cos  Or  cosh  0;  —  i  sin  0r  sinli  0/ )  .  (2  1 5) 

For  0/  negative,  Im[fcj]  is  positive,  and  e'klZ  represents  an  exponentially  decreas¬ 
ing  amplitude.  Introducing  the  change  of  variables  given  by  Eqs.  (2-13)  into 
Eq.  (2-12)  results  in 

ik  rf“*oo  r2ic 

p(x,  y,z)  =  —  /  /  sin  Od^dO  z>  0  ,  (2-16) 

I  IT  J$=Q  J<t>=0 

where  the  dependence  of  kx,  ky,  and  kz  on  0  and  <j>  is  summarized  in  Eqs.  (2-13). 

We  consider  now  a  point  source  S  at  coordinates  (0,0,  zs),  located  a  dis¬ 
tance  z„  above  a  plane  interface  separating  two  fluid  media  with  sound  speeds  c 
and  ci,  densities  p  and  pi,  and  wave  numbers  k  and  fcx,  respectively  (see  Fig.  2-2). 
The  field  at  a  given  receiver  R  at  horizontal  range  r  =  \/x7  -f  y2  and  height  zr 
above  the  interface  is  composed  of  a  direct  path  contribution  and  a  reflected  con¬ 
tribution.  The  reflected  pressure  field  p(r,zt)  at  R  is  constructed  by  multiplying 
each  of  the  plane  waves  in  the  integrand  of  Eq.  (2-16)  by  the  plane  wave  reflection 
coefficient  V (0)  and  by  taking  into  account  the  total  vertical  distance  zt  —  za  +  zr 
each  plane  wave  must  travel: 

P{r,  zt)  —  —  I  /  e,(fcxX+*»*+fc**‘)  V(0)  sin  Od<f>dO  .  (2  17) 

27T  Js= o  J<t>= a 

The  plane  wave  (Rayleigh)  reflection  coefficient  for  the  fluid-fluid  interface  is 


V(0)  = 


m  cos  0  -  \Zn2  —  sin2  0 
m  cos  0  -f  \/j)2  —  sin ^  0 


(2  18) 


where 


n  =  c/r,  ,  m  =  p\jp 


(2  19) 
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IMAGE 

SOURCE 


Figure  2-2  Geometry  FOR  THE  REFLECTED  FIELD. 

The  critical  angle  introduced  by  V(0)  is  defined  as 

0CT  =  arcsin(c/ci)  =  arcsin(n)  .  (2-20) 

The  integral  over  <}>  in  Eq.  (2-17)  can  be  performed  analytically  by  sub¬ 
stituting  for  kx  and  ky  from  Eqs.  (2-13)  and  expressing  the  receiver  coordinates 
x  and  y  in  terms  of  polar  coordinates  r  and  <)>i : 

x  =  rcos<f>\  y  =  rs\n<f>\  •  (2-21) 


The  result  is 


/*  =  /  ci(k'*+k*y)  d<j> 

J<t>= o 

,2* 

—  /  gikr  tin 0(co*4i  cos <t>+*\n <t>\ 

h=0 


_  f  eikr  sin  ))  ^ 

—  2nJ0(kr  sin  0)  . 

Introducing  this  result  in  Eq.  (2-17),  we  obtain 


(2-22) 


fi-IOO 

p(r,  zt)  =  ik  I  Jo{krs\n0)e,kziC<*eV(0)s\n0d0  .  (2-23) 

Jo 


12 


Equation  (2-23)  is  an  exact  representation  of  the  reflected  pressure  field. 
However,  in  order  to  apply  the  method  of  steepest  descent,  we  need  to  manipulate 
the  integral  once  more.  We  use  the  following  identity14  for  the  Bessel  function: 

Mu)  =  (l/2)lflJ"(u)  +  «S,,(«)l 

=  (1/2)  [/#>(..)- flf  >(-«)]  ,  (2-24) 

where  Hq  *  and  are  Hankel  functions  of  the  first  and  second  kind,  respec¬ 
tively,  of  order  zero.  After  substituting  Eq.  (2-24)  into  Eq.  (2-23),  we  obtain  two 
integrals  in  terms  of  Hq\  A  change  of  variable  from  6  to  ~0  is  then  introduced 
into  the  second  of  the  two  integrals.  Since  V(6)  and  cos&  are  even  functions  and 
sin#  is  an  odd  function,  the  resulting  integrand  is  the  same  as  that  of  the  first  in¬ 
tegral,  but  the  integration  runs  from  —n/2  +  ioo  to  0.  Therefore  the  two  integrals 
may  be  added  together  by  simply  combining  the  limits  of  integration: 

p(r,zt)  =  —  R{q  ] (kr  sin  0)eikztCoai  V{8)  sin  0dO  .  (2-25) 

2  J—  y  +  ioo 

The  contour  of  integration  in  the  complex  #-p!anc  is  shown  in  Fig.  2-3. 

The  first  approximation  introduced  into  Brekhovskikh’s  analysis  is  the 
use  of  the  asymptotic  expansion  for  the  Hankel  function 

H$\kr  sin  0)  ~  \  ■  eikr*mB  .  (2-26) 

V  ttAt  sin  0 

Substituting  for  the  Hankel  function  in  Eq.  (2-25),  we  obtain 

p(r,zt)  =  e^rsio$+z,coa9)V(0)y/^0d0  ,  (2-27) 

V  2irr  7-f+ioo 

valid  for  large  kr  sin  0.  Although  the  contour  of  integration  for  Eq.  (2-27)  passes 
through  0  =  0,  where  the  above  condition  clearly  is  not  satisfied,  we  shall  see  in 
Sec.  2-5  that  in  the  application  of  the  method  of  steepest  descent,  the  contour  of 
integration  is  deformed  in  such  a  way  that  it  does  not  pass  through  6  =  0.  The 
condition  of  large  fcrsin#  is  also  not  satisfied  at  short  ranges  r;  this  issue  will  be 
discussed  in  Sec.  2-7.2. 

Recognizing  from  Fig.  2-2  that  r  =  R\  sin#0  and  zt  =  Rx  cos  6o,  we  may 
write  the  reflected  pressure  field  as 

eikR,  co'(0-0o)  V(0)  y/^jj dd  ^ 


p{0o,R\)  = 


(2-28) 


CONTOUR  OF 
INTEGRATION 


BRANCH 


BRANCH 

POINTS 


Figure  2-3  Contour  of  integration  in  the  complex  0-plane  for  the  reflected 
field  integral  Eq.  (2-25).  Also  shown  are  the  branch  points  (at  ±0cr  =  ± arcsin(n)) 

AND  BRANCH  CUTS  INTRODUCED  BY  V{8). 
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where  the  receiver  position  is  now  specified  by  the  incident  angle  &q  and  the 
distance  R\  from  the  image  source  (see  Fig.  2-2).  Equation  (2-28),  finally,  is 
the  correct  form  of  the  integral  for  application  of  the  method  of  steepest  descent, 
which  will  be  described  next. 

2-3  THE  METHOD  OF  STEEPEST  DESCENT 
2-3,1  General  formulation 

The  method  of  steepest  descent  (see  Ref.  5,  Sec.  27)  is  used  to  approxi¬ 
mate  integrals  of  the  form 

I  =  f  epmF(0)d0  ,  (2-29) 

Jc 

where  C  is  a  general  contour  in  the  complex  0-plane,  p  (unrelated  to  the  density 
in  m)  is  real,  and  f(0)  and  F(0)  are  analytic  functions  of  the  complex  variable 
0.’  The  method  of  steepest  descent  consists  of  three  basic  steps:  (1)  finding  one 
or  more  saddle  points,  (2)  deforming  the  path  of  integration  into  the  steepest 
descent  path(s),  and  (3)  estimating  the  value  of  the  integral.  Saddle  points  are 
important  because  they  mark  particular  regions  of  the  complex  0-plane  where  the 
integrand  makes  a  distinct  contribution  to  the  integral.  The  path  of  integration  is 
deformed  into  the  path  of  steepest  descent  (which,  by  definition,  passes  through 
the  saddle  point)  because  the  contribution  to  the  integral  is  made  over  a  relatively 
short  interval.  To  obtain  an  estimate  of  the  integral,  the  integrand  is  replaced 
by  its  Taylor  series  expansion  about  the  saddle  point,  and  the  simpler  integral 
that  results  is  then  performed  analytically.  In  our  applications  of  the  method  of 
steepest  descent  to  acoustic  field  integrals,  we  find  that  saddle  points  represent 
eigenray  angles. 

Saddle  points,  which  we  will  call  7  throughout  this  work,  are  defined  by 
the  criterion 

f'{0)  =  0  at  0  =  7  ,  (2-30) 

where  prime  indicates  differentiation  with  respect  to  0.  Expressing  the  phase 
function  /  in  terms  of  its  real  and  imaginary  parts 

f(0)  =  M0)  +  if,(0)  ,  (2-31) 

we  define  the  steepest  descent  path  as  the  path  in  the  complex  0-plane  that  passes 
through  the  saddle  point  7.  and  along  which  //?  decreases  most  rapidly.  By  the 

"We  will  often  refer  to  f(9)  as  the  "phase  function,"  with  the  understanding  that  it  is  actually 
related  to  the  phase  of  the  exponential  by  the  constant  factor  i/p. 
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properties  of  the  real  and  imaginary  parts  of  an  analytic  function,  //  remains 
constant  along  the  path  of  steepest  descent  of  / r .  Since 


(2-32) 


we  see  that  along  the  steepest  descent  path  the  amplitude  of  the  exponential  de¬ 
creases  rapidly  as  one  moves  away  from  the  saddle  point,  while  the  phase  remains 
constant.  Mathematically,  the  steepest  descent  path  is  defined  by 

m  =  /(l)  -  ,  (2-33) 


where  $  is  a  real  number  that  parameterizes  the  path. 

To  obtain  an  approximation  to  the  integral,  we  equate  a  Taylor  series 
expansion  for  f(6)  about  the  saddle  point  7  to  the  steepest  descent  path  definition 
in  Eq.  (2-33): 

m  =  /(•»)- 

=  m  +  /'b)(»  - 1)  +  jA-rK*  - i)2  +  ■  •  •  .  (2  34) 

where  the  sec  nd  term  in  the  Taylor  series  is  zero  due  to  the  saddle  point  definition 
Eq.  (2-30).  Keeping  only  the  terms  shown  and  solving  for  s  in  terms  of  0,  we 
obtain  the  change  of  variable 


s  = 


(0-7) 


The  change  of  variable  in  Eq.  (2-29)  yields 


(2-35) 


'  (2'36) 

Note  that  the  integrand  is  small  at  all  values  of  s  except  in  the  neighborhood  of 
s  =  0.  If  we  assume  F(s)  is  slowly  varying  in  this  neighborhood,  F(s  =  0)  =  F( 7) 
can  be  brought  outside  the  integral,  leaving  only  the  definite  integral 

f  e-"1  ds  =  JTTp  .  (2-37) 

J -OO  V 


The  first  order  saddle  point  approximation  to  the  integral  is,  then, 

r  ■  ■  ■ 
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The  correct  sign  of  the  square  root  in  Eq.  (2-38)  is  found  by  considering  Eq.  (2-35) 
and  verifying  that  s  is  positive  in  the  forward  direction  of  integration  from  the 
saddle  point.  In  the  case  at  hand,  the  real  part  of  the  radical  should  be  positive. 
Brekhovskikh  derives  higher  order  approximations  in  Sec.  27  of  Ref.  5.  Additional 
terms  include  factors  of  p"^,  p-*,  and  so  forth. 


2-3.2  Classical  ray  theory 

A  straightforward  application  of  the  method  of  steepest  descent  to  the 
integral  for  the  reflected  field  in  Eq.  (2-28)  yields  a  result  consistent  with  classical 
ray  theory.  In  terms  of  the  parameters  given  in  Eq.  (2-29),  we  rewrite  Eq.  (2-28) 


p(0o,fli) 


where 


exp{(^Rij[i  cos(0  —  0O)] }  V (0)  '/sin  9  dO 

P  /(»)  F(6) 


(2-39) 


p  =  kR\  (2-40a) 

f(9)  =  icos(0-0o)  (2-40b) 

F{0)  =  y(0)VshT0  .  (2-40c) 


Upon  setting  /'{'))  =  — *sin(7  —  90)  =  0,  we  find  the  saddle  point  to  be  7  = 
0q.  Using  the  saddle  point  approximation  in  Eq.  (2-38)  and  recalling  that  r  = 
Ri  sin  0O,  we  write  the  reflected  pressure  field 


p{0o,Ri ) 


(2-41) 


which  represents  the  field  due  to  a  point  source  at  the  image  S'  (see  Fig.  2-2), 
whose  strength  is  multiplied  by  the  factor  V{0o)-  This  solution  expresses  the 
reflected  field  in  terms  of  spherically  spreading  rays  that  emanate  from  the  (ac¬ 
tual)  source  S  and  specularly  reflect  from  the  interface,  acquiring  the  plane  wave 
reflection  coefficient  factor  V(0O)  in  the  process. 

Comparisons  between  classical  ray  theory  [Eq.  (2-4  l)j  and  the  exact  so¬ 
lution  [obtained  by  numerically  integrating  Eq.  (2-23)]  at  various  frequencies  are 
shown  in  Fig.  2-4.  The  interface,  typical  of  a  water-sand  interface,  is  characterized 
by  n  =  0.866  (0CT  =  60.0°)  and  m  =  1.67.  The  vertical  axis  is  the  magnitude  of 
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Figure  2  4  Comparison  between  classical  ray  theory  [Eq.  (2-41)]  and  the  exact 
result  [Eq.  (2-23)]  for  the  normalized  reflected  field.  For  each  plot  of  the 
exact  field,  the  value  of  z,  (see  Fig.  2-2)  in  wavelengths  A  is  indicated.  Tiif. 
CRITICAL  ANGLE  LIES  AT  0cr  -  60.0°. 

the  reflected  pressure  field  in  decibels  (dB),  where  0  dB  corresponds  to  the  field 
at  unit  distance  from  a  point  source  in  free  space.  The  field  in  Fig.  2  1  has  been 
normalized  by  removing  the  spherical  wave  factor  c'kn' / R\  at  each  point.  All 
that  remains  is  the  effect  of  the  reflection  process,  which,  according  to  classical 
ray  theory,  is  just  V(0o)  [see  Eq.  (2  41)].  The  horizontal  axis  specifies  receiver 
location  by  giving  the  angle  of  incidence  0o  from  the  image  source,  where  :t  in 
wavelengths  X  is  constant  for  each  plot.  Note  that  t  in  classical  ray  result  does  not 
include  the  oscillations  in  the  exact  field  at  incident  angles  larger  than  critical. 
We  will  see  in  Sec.  2  5  that  these  oscillations  are  caused  by  interference  between 
the  totally  reflected  field  and  the  lateral  wave  field.  The  largest  error  in  classical 
ray  theory  occurs  at  the  lowest  frequency  (smallest  :,)  and  at  angles  (  lose  to  the 
critical  angle. 


d 


T 
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2-4  COMPLEX  RAY  FORMULATION 

2-4.1  Rearrangement  of  the  integral 

In  a  variation  on  the  method  of  steepest  descent  (Sec.  31.3,  Ref.  5), 
Brekhovskikh  allows  the  phase  of  the  reflection  coefficient  to  influence  the  locations 
of  the  saddle  points.  Expressing  V  in  terms  of  its  magnitude  IVj  and  phase  <f>v, 
we  may  write  Eq.  (2-28) 

p{0 o,  Ri )  =  i  fj-  e‘T  /7_,0°  ei[kR'  JV(6>)|  dO  .  (2-42) 

V  2 irr  J-f+.oo 

The  phase  function  in  the  exponent  of  Eq.  (2-29)  is  now 

f(Q)  =  i  cos (0  -  0O)  +  j^r>„(0)  >  (2-43) 

and  the  saddle  point  criterion  is 

f(0)  =  -i  sin(0  -  0O)  +  YrW)  =  0  at  0  =  7  .  (2-44) 

This  formulation  results  in  saddle  points  on  the  real  axis,  each  of  which  corre¬ 
sponds  to  the  angle  of  the  eigenray  that  starts  at  the  source  and  passes  through 
the  receiver.  It  is  the  formulation  used  in  Refs.  9-11. 

In  the  present  work,  we  take  Brekhovskikh’s  variation  one  step  further 
by  allowing  both  the  amplitude  and  phase  of  the  reflection  coefficient  to  affect  the 
saddle  point  location.*  To  do  this  we  express 

V{0)  ^e^v(9)  ,  (2-45) 


and  rewrite  Eq.  (2-28) 


p(OoJU) 


^kRi  [i  cos(0-0o)+  jTTj-  In  V (0)) 


\ZsirT?  d0 


(2-46) 


*A  careful  reading  of  Sec.  31.3  of  Ref.  5  suggests  that  Brekhovskikh  may  have  intended  just 
this  step.  He  sets  V(0)  —  c "*•(*)  and  states  that  “the  phase  coefficient  of  reflection  can  be 
a  complex  function.”  However,  he  docs  not  consider  complex  saddle  points,  and  later  in  the 
section  he  applies  this  formulation  only  to  the  case  of  total  reflection,  writing  the  phase  <f>v  as  a 
real  quantity. 


L 
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where  the  associations  with  Eq.  (2-29)  are  now 


p  =  kRi 

(2-47a) 

f(6)  =  icos(0  -  0O)  +  *^"ln  V(0) 

(2  -47b) 

F(9)  =  >/sin0  . 

(2-47c) 

The  saddle  point  criterion  is 

/'(*)= in(«-*„)  +  ji-^=0  =  7  ,  (2-48) 

and  the  second  derivative  /",  which  is  needed  for  the  saddle  point  approximation 
Eq.  (2-38),  is 

1  VI/" -H/')2 

f"{9)  =  ~tcos(e  -  0O)  +  — - •  (2-49) 

The  result  of  the  rearrangement  of  the  plane  wave  integral  and  reformulation  of 
the  saddle  point  criterion  is  that  the  saddle  points  7  are,  in  general,  located  at 
complex  values  of  0.  The  motivation  for  incorporating  V{0)  into  f(0)  is  that  we 
want  the  saddle  point  criterion  Eq.  (2-48)  to  include  the  complicated  effects  that 
the  reflection  process  introduces. 

The  saddle  point  approximation  to  the  integral  for  the  reflected  field  in 
Eq.  (2-46)  is  obtained  by  substituting  the  appropriate  quantities  into  Eq.  (2-38). 
After  simplification,  we  obtain 

*00,  a,)  =  vw  c'tR' '  (2-50) 

where  f"  is  given  in  Eq.  (2-49). 

2-4.2  Branch  cuts  in  V(0) 

The  dependence  of  the  plane  wave  reflection  coefficient  V  on  incident 
angle  0  is  what  makes  the  reflected  field  much  more  complicated  than  the  simple 
classical  ray  theory  result  in  Eq.  (2-41).  In  particular,  the  square  root  term 
q(0)  =  y/71?  —  sin2  0  =  \fq{0)  in  Eq.  (2-18)  for  V{0)  introduces  branch  points  and 
branch  cuts  in  the  complex  0-plane.  Branch  points  lie  at  values  of  0  for  which  the 
radicand  q  is  zero.  It  is  easy  to  see  that  the  branch  points  lie  at 


0  —  ±arcsin(u)  =  ±0cr  , 


(2  51) 


where  0CT  is  the  critical  angle.  The  branch  cuts  are  chosen  to  be  along  lines 
for  which  q(6)  lies  on  the  positive  real  axis.  This  imposes  the  following  two 
requirements: 

Re[g(0)]  =  njf  —  nj  —  (sin2  On  cosh2  0I  —  cos2  0fl  sinh2  0/)  >  0  (2-52a) 
Imf^#)]  =  2(n/jn/ +  sin  cos  0/i  sinh  0/ cosh  0/)  =  0  ,  (2-52b) 

where,  to  allow  for  the  possibility  of  attenuation  in  either  medium,  the  index  of 
refraction  is  written  as  the  complex  number  n  =  n#  +  inj.  For  the  case  of  no 
attenuation  (n/  =  0),  solution  of  Eqs.  (2-52a)  and  (2-52b)  for  0  results  in  the 
branch  cuts  shown  in  Fig.  2-3. 

Since  the  square  root  operator  is  double  valued,  the  evaluation  of  V{0) 
introduces  two  Riemann  sheets,  where 


lm[9(0)] 


>0  on  the  upper  Riemann  sheet 
<  0  on  the  lower  Riemann  sheet 


(2-53) 


The  branch  cuts  are  significant  in  our  problem  because  in  order  to  keep  the  func¬ 
tion  q(0 )  [and  therefore  V(0)]  continuous,  one  must  change  Riemann  sheets  when 
an  integration  path  crosses  a  branch  cut. 

The  contour  of  integration  for  the  reflected  field  integral  Eq.  (2-46)  begins 
on  the  upper  Riemann  sheet  at  0  =  — 7r/2  +  too  and  ends  on  the  upper  Riemann 
sheet  at  0  =  7r/2  —  ioo.  By  Cauchy’s  theorem  the  original  contour  may  be  de¬ 
formed  into  the  steepest  descent  path  as  long  as  no  singularities  are  crossed  in 
the  deformation  and  the  path  starts  and  ends  on  the  correct  sheets  at  the  correct 
points.  Therefore,  when  applying  the  method  of  steepest  descent,  it  is  essential  to 
verify  that  the  deformed  path  of  integration,  which  may  cross  the  branch  cuts  one 
or  more  times,  is  equivalent  to  the  original  contour.  When  the  steepest  descent 
path  of  a  saddle  point  ends  on  the  wrong  sheet,  we  must  locate  another  saddle 
point  whose  path  completes  the  required  contour. 


2-5  SADDLE  POINT  STRUCTURE  OF  THE  REFLECTED  FIELD 

The  saddle  points  and  steepest  descent  paths  for  the  reflected  field  along 
a  horizontal  line  are  shown  in  Fig.  2-5.  The  interface  is  characterized  by  0cr  =  60° 
and  m  =  1.67.  The  source- receiver  geometry  is  characterized  by  zt  =  5A,  and  0o  is 
varied  from  50°  to  80°  in  increments  of  2.5°.  The  saddle  point  angles  7  [obtained 
by  numerically  solving  Eq.  (2-48)]  as  a  function  of  receiver  position  (given  by  0O) 
are  indicated  by  the  bold  dots.  Each  saddle  point  represents  the  angle  at  which 
an  eigenray  arrives  at  the  receiver.  The  physical  picture  of  the  eigenray  paths  for 
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Figure  2-5  Saddle  points  and  steepest  descent  paths  for  the  reflected  field.  The 

GEOMETRY  IS  DESCRIBED  IN  THE  TEXT,  AND  THE  INTERFACE  IS  THE  SAME  AS  IN  FlG.  2-4 

Steepest  descent  paths  are  shown  for  the  saddle  points  in  the  region  of  the 

CRITICAL  ANGLE.  DASHED  PORTIONS  LIE  ON  THE  LOWER  RlEMANN  SHEET.  THE  INSET 
ILLUSTRATES  THE  COMPLETION  OF  THE  REQUIRED  CONTOUR  OF  INTEGRATION  FOR  THE  TWO 
BASIC  CASES. 


AS-89-231 
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both  the  reflected  and  transmitted  fields  is  shown  in  Fig.  2-6.  For  small  0q  (short 
ranges)  the  field  is  made  up  of  a  single  partially  reflected  eigenray.  For  large  0O 
(long  ranges)  the  field  is  made  up  of  a  totally  reflected  eigenray  plus  a  lateral  wave 
eigenray.  The  origin  of  the  lateral  wave  will  be  discussed  shortly. 


REFLECTED  FIELD 


Figure  2-6  Eigenray  picture  for  the  reflected  field  in  the  upper  medium  (see 
Sec.  2-5)  and  transmitted  field  in  the  lower  medium  (see  Sec.  2-8.2). 


Recall  that  classical  ray  theory  predicts  saddle  points  to  lie  at  7  =  0O. 
However,  from  Fig.  2-5  we  see  that  saddle  points  of  the  partially  reflected  eigenrays 
are  complex,  lying  slightly  below  and  to  the  left  of  0O,  and  that  saddle  points  of  the 
totally  reflected  eigenrays  lie  on  the  real  0-axis,  but  at  angles  noticeably  smaller 
than  0q.  In  Sec.  2-7  we  will  see  that  these  shifts  in  eigenray  angle  away  from  0O  are 
due  to  displacements  in  the  ray  trajectories  along  the  interface  upon  reflection, 
as  shown  in  Fig.  2-6.  Finally,  saddle  points  of  the  lateral  wave  eigenrays  lie  next 
to  the  critical  angle.  As  the  horizontal  range  to  the  receiver  increases,  the  lateral 
wave  saddle  points  move  closer  and  closer  to  the  critical  angle.  There  are  6  lateral 
wave  ray  saddle  points  in  Fig.  2-5  (to  match  the  6  totally  reflected  ray  saddle 
points),  but  since  they  all  lie  very  close  to  0cr,  it  is  not  possible  to  draw  them 
distinctly. 

Steepest  descent  paths  are  shown  for  the  saddle  points  in  the  region  of 
the  critical  angle.’  The  paths  of  the  first  seven  saddle  points,  corresponding  to 


The  paths  drawn  in  Fig.  2  -5  are  those  along  which  the  numerical  integration  technique, 


partially  reflected  eigenrays,  cross  the  branch  cut  on  the  real  0-axis.  Since  these 
paths  cross  the  branch  cut  again  on  the  positive  imaginary  axis  (see  inset  of 
Fig.  2-5),  they  begin  on  the  correct  (upper)  Riemann  sheet  and  thereby  complete 
the  required  contour.  However,  the  path  of  the  first  totally  reflected  eigenray, 
for  the  receiver  at  0O  =  67.5°,  does  not  cross  the  branch  cut  on  the  real  axis,  so 
it  begins  on  the  wrong  (lower)  Riemann  sheet.*  At  this  point  it  is  necessary  to 
seek  an  additional  saddle  point  (an  additional  root  of  Eq.  (2-48)]  whose  steepest 
descent  path  completes  the  contour.  Such  a  saddle  point  is  found  on  the  real  0-axis 
close  to  the  branch  point  0cr  =  60°.  The  saddle  point  represents  the  contribution 
from  the  lateral  wave,  and  its  steepest  descent  path  loops  around  the  branch  point 
such  that  it  starts  on  the  upper  Riemann  sheet  and  ends  on  the  lower  Riemann 
sheet  (see  inset  of  Fig.  2-5).  Thus,  when  the  lateral  wave  ray  and  totally  reflected 
ray  are  added  together,  their  integration  paths  complete  the  required  contour. 

The  physical  interpretation  of  the  lateral  wave,  illustrated  in  Fig.  2-6, 
is  well  known  (see,  for  example,  Sec.  30.2  of  Ref.  5).  It  is  generated  by  energy 
incident  on  the  interface  in  the  region  of  the  critical  angle.  Some  of  this  energy 
penetrates  into  the  bottom  and  travels  almost  parallel  to  the  interface  at  the 
bottom  sound  speed.  Since  the  pressure  across  the  interface  must  be  continuous, 
the  horizontally  traveling  energy  in  the  bottom  continuously  excites  a  wave,  the 
lateral  wave,  that  travels  at  the  critical  angle  in  the  water.  The  lateral  wave 
arrives  at  all  receivers  that  lie  beyond  a  certain  range  from  the  source,  and  its 
wavefronts  are  conical  about  the  S-S'  axis  in  Fig.  2-6. 

The  transition  from  a  field  composed  of  one  saddle  point  at  a  complex 
angle  (representing  a  partially  reflected  eigenray)  to  a  field  composed  of  two  saddle 
points  at  real  angles  (representing  the  totally  reflected  eigenray  and  the  lateral 
wave  eigenray)  occurs  at  a  ray  theory  caustic.  Caustics,  which  can  be  caused  by 
the  bending  of  rays  as  they  are  continuously  refracted  by  a  sound  speed  gradient 
in  the  medium,  are  points  in  space  where  two  rays  coalesce  into  one.  In  the  present 
case,  the  caustic  is  caused  by  the  reflection  process,  where  two  separate  rays  are 
formed  upon  reflection.  When  only  considering  rays  with  real  angles,  the  caustic 
is  the  transition  point  between  the  shadow  zone,  into  which  no  rays  penetrate, 
and  the  insonified  zone,  where  two  rays  arrive  at  every  point.  Mathematically, 
caustics  are  defined  as  saddle  points  [/'(-y)  =  0]  at  which  /"(y)  =  0. 

In  Fig.  2-5  the  caustic  saddle  point  is  recognized  as  the  point  where  the 
totally  reflected  rays  merge  with  the  lateral  wave  rays  (7,.  %  62°).  The  location 

described  in  Sec.  2-G,  is  performed. 

*  Although  the  transition  from  total  to  partial  reflection  is,  in  reality,  a  gradual  process,  we 
use  this  criterion  of  steepest  descent  paths  crossing,  or  not  crossing,  the  branch  cut  to  distinguish 
between  partially  and  totally  reflected  rays 
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of  the  caustic  in  physical  space  is  given  by  zt  =  5A  and  Oo  ss  66.5°  (obtained  by 
interpolating  between  the  value  of  &o  =  65°  for  the  last  partially  reflected  ray 
shown  and  0O  =  67.5°  for  the  first  totally  reflected  ray  shown).  Since  /"(7c)  =  0, 
the  saddle  point  approximation  for  the  field  given  by  Eq.  (2-50)  is  infinite  at  the 
caustic.  In  the  following  subsection  we  summarize  our  method  for  finding  the  field 
in  such  regions  where  Eq.  (2-50)  is  invalid. 

Attenuation  in  the  bottom  can  be  included  by  introducing  an  imaginary 
part  to  the  sound  speed  Cp*  As  a  result,  the  index  of  refraction  n  =  c/c,  in 
Eq.  (2-18)  is  complex,  the  critical  angle  0CT  =  arcsin(n)  acquires  a  positive  imag¬ 
inary  part,  and  the  saddle  points  are  shifted.  Figure  2-7  is  the  same  as  Fig.  2-5, 
except  the  lower  medium  has  an  attenuation  of  oaB^O.S  dB/A.  Notice  that  in  this 
case  there  is  a  smooth  transition  between  the  partially  reflected  and  the  totally 
reflected  ray  saddle  points.  Since  there  is  no  receiver  location  where  two  saddle 
points  coalesce  into  one,  there  is  no  caustic.  The  lateral  wave  saddle  points  have 
followed  Oct  (the  branch  point)  off  the  real  axis,  and  the  corresponding  eigenrays 
suffer  attenuation  as  they  travel  along  the  interface  in  the  bottom.  The  totally 
reflected  ray  saddle  points  have  a  small  negative  imaginary  part,  and  the  eigenrays 
are  attenuated  to  a  lesser  degree. 

2-6  COMPUTATION  OF  THE  EIGENRAY  FIELD 

The  saddle  point  approximation  in  Eq.  (2-50)  is  based  on  the  assumption 
that,  in  the  neighborhood  of  7,  the  phase  function  f(0)  is  quadratic  and  the 
steepest  descent  path  is  a  straight  line  in  the  complex  0-plane.  But  this  assumption 
is  invalid  when  there  is  significant  curvature  in  the  steepest  descent  path  (see,  for 
example,  Fig.  2-5).  At  the  caustic  itself,  Eq.  (2-50)  predicts  an  infinite  eigenray 
magnitude. 

One  way  to  find  the  correct  field  in  the  caustic  region  would  be  to  apply 
caustic  corrections  that  are  based  on  uniform  asymptotic  expansions.11,16  However 
this  approach  is  not  attractive  for  several  reasons.  First,  it  is  difficult  to  smoothly 
join  the  shadow-zone  field  with  the  field  due  to  the  partially  reflected  rays.  Second, 
the  asymptotic  expansions  become  extremely  cumbersome  as  one  is  forced  to 
include  higher  order  terms.  In  our  case,  low  frequencies  and  the  proximity  of  the 
caustic  to  the  singularity  at  the  critical  angle  would  make  the  calculation  of  higher 
order  terms  necessary. 

‘This  procedure  is  summarized  in  Sec.  A8.4  of  Ref.  15.  We  define  a  complex  sound  speed 
c  =  c/i  +  ic/.  For  small  attenuations,  we  set  the  imaginary  part  cj  =  —[(In  10)/(40T)]cfla<jB, 
where  cr  is  the  usual  (real)  sound  speed  and  q^b  is  the  attenuation  in  the  medium  in  units 
db/A. 


Figure  2-7  As  in  Fig.  2-5,  except  n  =  0.866  +  »0. 009,  corresponding  to  an  attenua¬ 
tion  of  0.5  dB/A  in  the  lower  medium. 
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To  find  the  eigenray  field  for  saddle  points  in  the  region  of  the  critical 
angle  we  simply  integrate  numerically  along  the  steepest  descent  path.  The  ad¬ 
vantages  of  this  technique  are  that  it  works  at  any  frequency,  it  can  be  used  to 
determine  when  more  than  one  saddle  point  must  be  found,  and  it  can  be  applied 
to  any  integrand  for  which  the  phase  function  /  and  its  derivatives  /'  and  /"  can 
be  computed. 

Points  along  the  steepest  descent  path  are  found  sequentially,  working 
away  from  the  saddle  point  in  one  direction,  then  the  other.  The  procedure  is  to 
determine  desired  values  of  /  along  the  path,  and  then  use  a  Taylor  series  for  / 
to  obtain  an  initial  guess  for  each  successive  6  point  on  the  path.  The  first-order 
Taylor  series  at  the  previous  point  6P  on  the  steepest  descent  path  is 

m  =  m) + noP)(0  -  ep)  .  (2-5i ) 

Given  a  desired  value  of  f(9)  at  the  next  point  on  the  path,  fd ,  the  initial  guess, 
9g ,  is  obtained  by  solving  Eq.  (2-54): 

».-»,  + Ml -/(M/m)  •  (2-55) 

At  the  saddle  point  itself,  0P  =  7  and  f(0p )  =  0,  so  a  second-order  Taylor  expan¬ 
sion  is  needed.  The  result  is 

0,  =  7  +  \j‘l[U-f{'l)}ini)  ,  (2-56) 

where  the  real  part  of  the  square  root  term  is  chosen  to  be  positive  in  the  forward 
direction  of  integration  and  negative  in  the  opposite  direction. 

The  remaining  question  is  how  to  choose  the  desired  values  /j.  The 
imaginary  part  of  fd  must  be  held  constant  at  //(7),  the  value  at  the  saddle  point, 
in  order  to  keep  the  phase  of  the  exponential  constant  at  pf/( 7)  [see  Eq.  (2-32)]. 
The  real  part  of  fd  is  chosen  so  as  to  sample  the  amplitude  of  the  exponential  factor 
at  desired  levels,  referenced  to  the  maximum  of  i4mox  =  at  the  saddle  point. 

For  example,  if  the  first  sample  away  from  the  saddle  point  is  desired  to  have  an 
amplitude  of  0.99Amai,  then  we  would  set 

Mfd\  =  /r(  7)  +  ln(  .99)/ p  .  (2-57) 

The  desired  levels  at  which  to  sample  the  amplitude  factor  are  chosen  so  that  a 
linear  interpolation  between  samples  gives  a  good  approximation  to  a  Gaussian, 
the  approximate  shape  of  the  amplitude  function  along  the  steepest  descent  path. 
We  use  a  total  of  37  sample  points  in  the  numerical  integration,  resulting  in  an 
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error  in  the  pressure  field  on  the  order  of  ±0.01  dB.  The  end  points  of  the  inte¬ 
gration  paths  drawn  in  Figs.  2-5  and  2-7  are  determined  by  the  criterion  that  the 
exponential  factor  have  an  amplitude  of  0.005/4maI  there.  The  integration  algo¬ 
rithm  must  include  special  provisions  for  sections  of  paths  that  change  direction 
rapidly,  such  as  those  associated  with  the  lateral  wave,  which  loop  tightly  around 
the  critical  angle. 


The  results  using  numerical  integration  along  the  steepest  descent  path 
for  the  cases  in  Fig.  2-4  are  indistinguishable  from  the  exact  results  shown.  For  the 
case  of  Zt  —  5A,  the  steepest  descent  method  is  approximately  twice  as  efficient  as 
the  direct  numerical  integration  of  Eq.  (2-23).  However,  the  latter  requires  more 
sample  points  as  the  frequency  is  increased,  while  the  number  of  sample  points 
for  the  steepest  descent  method  is  independent  of  frequency.  When  the  saddle 
point  approximation  can  be  used  instead  of  steepest  descent  path  integration,  the 
efficiency  of  computing  the  field  improves  greatly. 


It  is  difficult  to  determine  from  the  saddle  point  characteristics  alone  how 
accurate  the  saddle  point  approximation  is  in  a  particular  case.  The  method  we 
employ  is  to  begin  in  the  area  where  the  approximation  is  the  poorest,  the  region 
of  the  critical  angle,  and  work  away  from  it  in  both  directions.  At  each  point 
we  compare  the  result  from  steepest-descent  integration  with  the  saddle  point 
approximation.  When  the  latter  is  sufficiently  close  to  the  “exact”  answer,  we 
do  not  need  to  perform  the  numerical  integration  for  any  more  receivers  in  that 
direction. 


The  accuracy  of  the  method  of  numerically  integrating  along  the  path  of 
steepest  descent  begins  to  deteriorate  only  when  the  frequency  becomes  extremely 
low.  For  zt  =  1A,  the  steepest  descent  method  differs  from  the  exact  solution 
[obtained  by  numerical  integation  of  Eq.  (2-23)]  by  less  than  0.1  dB  at  all  angles. 
For  zt  <  0.5A,  we  find  that,  instead  of  proceeding  directly  to  the  end  point  of  the 
integration  path  at  tt/2  —  too,  the  steepest  descent  paths  can  end  at  the  angle  of 
intromission,  6-,m,  where  K(0jm)  =  0  and  Re[/]  approaches  — oo.  For  the  interface 
considered  in  Fig.  2-5,  0,m  =  90°  —  i21°.  In  such  cases,  another  saddle  point  that 
lies  near  0jm  must  be  found.  Its  steepest  descent  path  begins  at  6im  and  ends  in  the 
correct  direction,  thus  completing  the  contour  of  integration.  In  Sec.  2-8.2  and 
Fig.  2-12  we  encounter  a  similar  phenomenon  for  the  transmission  case,  caused 
by  the  fact  that  the  transmission  coefficient  is  zero  at  a  particular  angle. 


J 
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2-7  GEOMETRIC  INTERPRETATION 

2-7.1  Derivation  of  eigenray  characteristics 

One  of  the  most  attractive  features  of  the  saddle  point  approach  is  that 
the  mathematical  solution  can  be  interpreted  in  terms  of  intuitive  ray  quantities, 
such  as  the  path  of  propagation,  the  geometric  spreading  loss,  the  reflection  loss, 
and  the  travel  time.  By  moving  the  reflection  coefficient  V(0)  into  the  phase 
function  /  in  Eq.  (2-46),  we  have  incorporated  the  reflection  process  into  the 
calculation  of  these  ray  characteristics. 

The  geometric  interpretation  of  the  ray  path  is  shown  in  Fig.  2-8.  The 


Figure  2-8  Eigenray  interpretation  of  the  saddle  point  7  for  the  reflected 
FIELD 

eigenray  at  angle  7  is  traced  from  the  source  S  to  the  interface  at  A  and  backward 
from  the  receiver  R  to  the  interface  at  B.  The  resulting  ray  displacement  AB  is 
defined  as  A.  Although  the  angle  7  is,  in  general,  complex,  we  can  only  draw 
it  as  a  real  angle.  In  cases  where  7  is  complex,  we  must  use  trigonometric  func¬ 
tions  defined  for  complex  arguments  in  the  derivations  that  follow.  The  physical 
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interpretation  of  eigenrays  with  complex  angles  will  be  discussed  in  Sec.  2-7 .3. 

We  begin  by  deriving  the  general  form  for  the  ray  displacement  A.  Using 
the  geometrical  relation 

sin(7  —  )  =  -S'D/Ri  =  —  Acos7//?i  (2-58) 


in  the  saddle  point  criterion 


1  V'iVi 

/'( 7)  =  -*  sin(7  -  0o)  +  Tn-  =  0 


we  obtain 


A(7)  = 


kRx  Vin) 

v\l) 


(2-59) 


(2-60) 


kc os  7  ^(7) 

Equation  (2-60)  reduces  to  Brekhovskikh’s  result  in  Eq.  31.30  of  Ref.  5  for  the 
special  case  where  7  is  real  and  | Vj  =  1  (total  reflection).  In  this  case,  A  = 
— <f>'v/(k  cos*/),  where  <f>v  is  the  phase  of  V. 

We  can  convert  the  saddle  point  criterion  Eq.  (2-59)  to  an  eigenray  equa¬ 
tion  by  (1)  substituting  for  V'/V  in  terms  of  A,  (2)  expanding  the  sin(7  —  0o) 
term,  and  (3)  using  the  geometric  relations  R\  sin0o  =  r  and  Rx  cos0o  =  z,  +  zr: 


/'( 7)  =  —  i(sin7cos0o  —  sin  0O  cos 7)  +  — jr-  (—  ikAcos^)  —  0 


(— )  /'( 

\c0s7  ) 


kRx 


'(7)  =  Rx  cos  0Q  tan  7  -  R\  sin  0O  +  A 

0  =  {z,  tan  7  +  A  -f  zr  tan  7)  —  r 
By  Fig.  2-8,  we  interpret  Eq.  (2-62)  as  the  eigenray  equation 

r  =  r,(7)  + A(7)  +  r2(7)  . 


(2-61) 

(2-62) 

(2-63) 


Equation  (2-63)  gives  the  horizontal  range  r  covered  by  a  ray  launched  at  complex 
angle  7.  Note  that  any  of  the  three  intermediate  range  terms  r1(  A,  and  r2  shown 
in  Fig.  2-8  may  be  complex,  but  their  sum  must  be  the  real  quantity  r. 

Upon  using  the  geometric  relation  (see  Fig.  2-8) 

Rx  cos(7  —  0O)  =  DC  +  CB  +  DR  =  ISA  -f  A  sin  7  +  BR  (2-64) 


for  the  phase  factor  in  the  saddle  point  approximation  Eq.  (2-50),  we  obtain  the 
reflected  field  in  terms  of  eigenray  characteristics: 


p(*o,/?i) 


exp(z  k(§~A  +  A  sin  7  +  BR)} 
ray  path  phase  $ 


(2-65) 
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The  factors  G+  and  G$  are  geometric  spreading  factors  and  will  be  described  in 
detail  in  the  next  section.  The  reflection  loss  factor  in  Eq.  (2-65)  is  just  the  plane 
wave  coefficient.  The  ray  path  phase  4>  shows  that  the  ray  travels  along  segments 
5 A  and  TTR  at  the  water  sound  speed  c,  but  along  the  interface  at  the  speed 

c&  =  c/  sin  7  .  (2-66) 

Thus,  lateral  wave  rays  (see  Sec.  2-5),  for  which  sin  7  «  sin0cr  =  c/cj,  travel  the 
distance  A  along  the  interface  at  the  bottom  sound  speed,  while  rays  at  grazing 
incidence,  7  =  90°,  travel  at  the  water  sound  speed.  Totally  reflected  rays,  for 
which  0CI  <  7  <  90°,  travel  at  speeds  somewhere  between  c  and 

2-7.2  Analysis  of  the  geometric  spreading  factors 

The  first  two  factors  in  Eq.  (2-65),  and  Go,  are  the  geometric  spread¬ 
ing  factors  in  the  two  orthogonal  directions  <f>  and  9  shown  in  Fig.  2-1.  The 
spreading  of  a  ray  bundle  in  the  0-direction  is  illustrated  in  Fig.  2-9(a),  which 
is  drawn  in  the  plane  of  incidence  of  the  ray  (the  plane  containing  the  ray  and 
the  normal  to  the  interface).  The  bundle  has  width  A 0  and  is  centered  about 
the  eigenray  angle  7.  Its  cross  section  at  the  receiver  is  r'(7)  cos  7A0,  where  r' 
is  the  derivative  of  r  in  Eq.  (2-63)  with  respect  to  7.  Since  the  geometric  loss  is 
inversely  proportional  to  the  square  root  of  the  variation  in  the  cross  section  of 
the  bundle,  we  expect  Go  —  ^/l/r'(7)  cos  7.  This  is  indeed  equivalent  to  the  ex¬ 
pression  in  Eq.  (2-65)  because  differentiation  of  Eq.  (2-61)  results  in  the  relation 
iRif"( 7)  =  r'(7)cos7.  Carrying  out  the  differentiation,  one  can  express 


r'( 7)  cos  7 

where 

A,  i  vv"-{vy 

A'  =  7 - TJY  ■  +  A  tan  7  .  (2-68) 

k  cos  7  Vi 

Near  the  caustic  Eq.  (2-67)  is  invalid  because  r'  approaches  zero,  and  G$  therefore 
approaches  infinity.  In  this  case,  the  exact  field  p  must  be  obtained  by  numerical 
integration  .dong  the  steepest  descent  path.  If  desired,  the  value  of  Ge  may  be 
obtained  by  rearranging  Eq.  (2-65): 


SA  +  A'  cos  7  +  BR 


(2-67) 


p(9q,  R\  ) 

(<4,  V{~t)  exp(<<t>) 


(2-69) 
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The  spreading  of  a  ray  bundle  in  the  ^-direction  is  illustrated  in  Fig. 
2-9(b).  A  ray  bundle  of  angular  width  A<£,  when  projected  on  the  plane  of  the 
interface,  has  angular  width  A<f>/  sin  7.  Cylindrical  symmetry  dictates  that  the 
center  of  curvature  of  the  various  portions  of  the  ray  path  lie  on  the  line  through 
the  source  5  and  perpendicular  to  the  interface.  Thus,  the  (linear)  cross  section 
of  the  bundle  at  the  receiver  is  simply  rA<£/  sin  7.  In  terms  of  the  path  lengths 
along  each  segment  of  the  ray  trajectory,  can  be  expressed  as 


— 


/sin  7 


SA  +  A/  sin7  +  BR 


(2-70) 


Note  that  when  the  reflection  coefficient  is  not  a  function  of  incident  angle,  A  =  0 
and  Go  =  G $  =  \J \  j R\.  In  such  a  case,  the  classical  ray  result  Eq.  (2-41)  is  exact. 

At  this  point  it  is  worthwhile  to  recognize  the  origin  of  the  spreading  loss 
factor  G#.  The  spreading  in  the  (^-direction  can  be  traced  back  to  the  change 
of  variable  in  Eq.  (2-14).  The  first  factor  on  the  right  hand  side  of  Eq.  (2-14) 
accounted  for  the  ^-spreading  and  introduced  the  sin  9  factor  inside  the  double 
integral  of  Eq.  (2-16).  The  integral  over  (f>  eventually  resulted  in  the  Hankcl 
function  factor  in  Eq.  (2-25),  which  was  then  replaced  by  its  asymptotic  form  in 
Eq.  (2-26).  Combination  of  the  asymptotic  form  of  the  Hankel  function  with  the 
sin  0  factor  resulted  in  the  factor  ^/sin  0/r. 

When  viewed  geometrically  as  in  Fig.  2-9,  it  makes  sense  to  separate  the 
spreading  factors  in  the  two  independent  directions.  This  is  achieved  by  evaluating 
the  ^-spreading  factor  G4,  =  \Js'\n9/r  at  the  saddle  point  7  and  taking  it  outside 

the  integral  over  6  in  Eq.  (2-46),  which  accounts  for  the  ^-spreading  factor  Go- 
With  the  Gj,  factor  removed,  the  Hankel  function  approximation  ( kr  sin  0  not 
small)  no  longer  causes  errors  at  short  ranges  and  low  frequencies.  As  we  shall 
see  in  Chapter  3,  separation  of  the  two  spreading  factors  in  this  way  enables 
our  basic  method  for  evaluating  the  field  to  be  extended  to  more  complicated 
environments,  where  ray  paths  may  have  multiple  reflections  at  possibly  different 
angles  of  incidence. 


2-7.3  Interpretation  of  complex  eigenrays 

As  we  have  seen,  saddle  points,  and  therefore  eigenray  angles,  can  be 
complex.  When  7  is  complex,  the  eigenray  characteristics  in  Eq.  (2-65)  have  non¬ 
zero  imaginary  parts  that  may  not  have  a  clear  physical  interpretation.  When 
saddle  points  are  complex  due  to  attenuation  in  one  or  both  media,  as  is  the 
case  for  the  saddle  points  representing  the  totally  reflected  and  lateral  wave  rays 
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in  Fig.  2-7,  the  imaginary  parts  of  the  eigenray  characteristics  may  be  viewed 
as  simply  accounting  for  the  attenuation  along  the  eigenray  paths.  However, 
even  with  no  attenuation,  partially  reflected  eigenray  angles  have  small  imaginary 
parts,  as  seen  in  Fig.  2-5.  The  largest  imaginary  parts  of  7  occur  at  receivers  near 
Oq  =  0CT.  The  significance  of  the  complex  saddle  points  in  this  region  is  that  the 
field  does  not  behave  like  a  simple  (classical  ray  theory)  ray  field  there.  At  lower 
frequencies  (source  and  receiver  closer  to  the  interface  in  terms  of  wavelengths)  the 
field  in  this  region  behaves  less  like  a  simple  ray  field,  and  the  imaginary  parts  of 
the  saddle  points  7  are  correspondingly  larger.  Although  we  can  assign  no  simple 
physical  interpretation  to  comp  lex- valued  eigenray  characteristics,  we  can  deduce 
the  approximate  path  of  the  energy  the  “ray”  represents  from  the  real  part  of  the 
eigenray  angle. 

The  possibility  of  complex- valued  saddle  points  makes  the  process  of  “ray 
tracing”  more  difficult.  An  eigenray  with  a  complex  angle  begins  and  ends  in 
real  space  at  the  source  and  receiver,  respectively,  but  otherwise  travels  through 
complex-valued  spatial  coordinates.  Thus,  in  regions  where  eigenrays  are  complex, 
we  may  not  use  the  usual  practice  of  launching  a  series  of  real  rays  and  tracing 
them  through  real  space. 

2-8  THE  TRANSMITTED  FIELD  IN  THE  BOTTOM 

In  this  section  we  solve  for  the  transmitted  field  in  the  bottom  using  the 
same  method  used  for  the  reflected  field.  After  developing  the  plane  wave  integral 
for  the  field  in  the  bottom,  we  summarize  the  saddle  point  structure  and  interpret 
the  field  in  terms  of  eigenray  characteristics. 

2-8.1  Integral  representation  of  the  transmitted  field 

The  derivation  of  the  plane  wave  integral  for  the  field  at  a  point  in  the 
bottom  begins  with  the  representation  of  the  spherical  wave  in  Eq.  (2-16).  A 
plane  wave  incident  on  the  interface  at  angle  9  is  transmitted  into  the  bottom  at 
angle  9X  according  to  Snell’s  law 

sin#i=sin0/n  (2-71) 

Figure  2-14  shows  the  geometrical  configuration  of  the  source  and  receiver.  We 
express  distances  in  the  bottom  as  well  as  in  the  water  in  terms  of  the  wavelength 
A  in  the  water. 

The  pressure  field  in  the  bottom  is  constructed  by  multiplying  each  plane 
wave  in  the  integrand  of  Eq.  (2-16)  by  the  transmission  coefficient  W(0)  and  by 
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SOURCE 


Figure  2-10  Geometry  for  calculation  of  the  transmitted  field  in  the  bottom 


taking  into  account  the  total  phase  in  the  water  kzza  and  in  the  bottom  kizzr : 

p(r ,  zr)  =  —  /T_,0°  e'(k*x+k»v+k‘Zt+ki’*r)  W(9)  sin  0  d<f>  dO  ,  (2-72) 

2ir  Je= o  J<t>= o 


where  kx,  ky,  and  kz  are  given  in  Eqs.  (2-13),  k}z  =  kx  cos#!  is  the  vertical 
component  of  the  wave  number  in  the  bottom,  and  the  fluid -fluid  plane  wave 
transmission  coefficient  is 


W{0)  = 


2m  cos  9 

m  cos  9  +  y/n 2  —  sin2  9 


(2-73) 


Following  the  derivation  for  the  reflection  case,  we  perform  the  integration  over  <j> 
analytically,  with  the  result 


p(r,z„zT)  =  ik  r  J0(kr  sin  9)  e'(k‘z‘+kl‘Zr)  W{9)  sin  9  d9  .  (2-74) 

Jo 


Next,  we  express  the  Bessel  function  in  terms  of  Hankel  functions  and  use  the 
asymptotic  form  of  the  Hankel  function.  Finally,  after  rearrangement  such  that 
W{9)  appears  in  the  phase  function,  we  obtain 


p{90,Rl,zr) 


\/-^-e'T  f 2  exp  {fcR,  (icos(0  -  0O) 
V  27 rr  J_*+I<x> 

-f  y-r-  In  W{9)  \  }  \/sin  9 d9 

k  it  | 


—  (cos  6  —  n  cos  i ) 

R\ 


(2-75) 


35 


The  integration  contour  and  branch  cuts  are  the  same  as  for  the  reflection  case 
in  Fig.  2-3.  The  choice  of  the  Riemann  sheet  [see  Eq.  (2  53)]  applies  to  the  term 
n  cos  =  Vn2  —  sin^  0  =  q  in  the  exponent  of  Eq.  (2  -75)  and  in  W(0).  Note  that 
the  transmitted  field  depends  on  zr  as  well  as  90  and  R\ . 

The  phase  function  in  the  transmitted  field  integral  Eq.  (2-75)  is 

f(0)  =  icos(0  -do)  —  l-^-(cos6  -  ncos0i)  +  — —  lnW(0)  ,  (2-76) 

kill 

the  saddle  point  criterion  is 

/'(0)  =  -isin(0-0o)  +  ^:(sin0  +  g/)  +  ^-^^  =  O  at  0  =  7,  (2-77) 

and  the  second  derivative  required  by  the  saddle  point  approximation  is 

x  izT  „  1  WW"  -  (W)2 

5  (0)  =  -icos(0 -0O)  +  —  {cos0  +  q  )  +  — - -  ,  (2-78) 

where 

q  =  —  sin0cos0/(ncos0i) 

q"  =  [sin2  0  —  cos2  0  —  (q')2]/(n  cos0i)  . 

2-8.2  Saddle  point  structure  of  the  transmitted  field 

The  saddle  points  and  steepest  descent  paths  are  found  in  a  manner 
analogous  to  the  reflection  case.  In  Fig.  2-11  we  show  the  situation  in  the  complex 
0-plane  for  a  low  frequency  case  where  za  =  4A,  zT  —  1  A,  and  the  angle  0O  specifying 
the  receiver  position  is  varied  from  50°  to  80°  in  increments  of  2.5°.  We  see 
that  for  the  seven  receivers  closest  to  the  source,  the  steepest  descent  path  for 
the  transmitted  ray  completes  the  required  contour  of  integration  by  itself.  But 
for  receivers  farther  from  the  source,  the  path  turns  back  up  and  to  the  left 
instead  of  down  and  to  the  right.  For  these  receivers  the  steepest  descent  path  of 
an  additional  saddle  point,  corresponding  to  the  evanescent  field,  completes  the 
contour.  The  evanescent  field  is  produced  by  energy  that  is  totally  reflected  at  the 
interface,  but  that  nevertheless  causes  a  disturbance  in  the  bottom.  The  complex 
angle  of  the  “evanescent  ray”  has  a  real  part  that  is  slightly  greater  than  0O  and 
an  imaginary  part  that  is  positive.  The  eigenray  interpretation  of  the  field  in  the 
bottom  is  shown  in  Fig.  2  (>. 

When  the  source  and  receiver  are  very  close  to  the  interface  and  some¬ 
what  far  apart,  a  complication  in  the  steepest  descent  path  of  the  evanescent  ray 
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6  R-  deg 


Figure  2-11  Saddle  points  and  steepest  descent  paths  for  the  transmitted  field. 
The  geometry  is  described  in  the  text,  and  the  interface  is  the  same  as  in  Fig.  2-4 
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saddle  point  can  arise.  Instead  of  ending  in  the  correct  direction,  the  path  ends  at 
0  =  90°  +t0°,  where  W  by  Eq.  (2-73)  is  zero  and  Re[/]  by  Eq.  (2-76)  approaches 
— oo.  When  this  happens,  we  must  search  for  an  additional  saddle  point  7eu2  whose 
path  of  integration  completes  the  required  contour  of  integration  (see  Fig.  2-12). 
We  call  7„.?,  which  lies  slightly  below  and  to  the  right  of  90°  in  the  complex 


Figure  2-12  SADDLE  POINTS  AND  STEEPEST  DESCENT  PATHS  FOR  A  CASE  WHERE  AN  AD¬ 
DITIONAL  “EVANESCENT  RAY”  SADDLE  POINT  MUST  BE  FOUND.  The  PARAMETERS  ARE: 
n  =  0  866,  m  =  1.67,  z,  -  zr  -  0.25A,  r  =  40A. 

0-plane,  a  second  “evanescent  ray”  saddle  point  because  it  represents  an  inhomo¬ 
geneous  plane  wave  that  travels  parallel  to  the  interface  (Re[7ew2]  ~  90°)  and  is 
exponentially  decreasing  in  the  bottom  (Im[7e„2J  <  0).  For  the  case  in  Fig.  2  -12, 
the  field  due  to  7e„2  is  just  5  dB  weaker  than  that  due  to  the  usual  evanescent 
field  saddle  point,  so  the  contribution  of  the  former  cannot  be  ignored.  We  do  not 
believe,  however,  that  the  additional  saddle  point  represents  any  distinct  mech¬ 
anism  of  propagation,  but  simply  that  it  is  sometimes  required  to  complete  the 
description  of  the  evanescent  field. 

It  is  interesting  to  note  how  the  saddle  point  structures  of  the  reflected 
and  transmitted  fields  are  related  at  the  interface,  where  the  total  pressure  field 
must  be  continuous.  At  short  ranges  (small  0O),  the  saddle  point  corresponding 
to  the  partially  reflected  ray  becomes,  as  the  interface  is  crossed,  the  saddle  point 
corresponding  to  the  transmitted  ray.  At  long  ranges  (large  0O),  it  is  the  reflected 
lateral  wave  ray  that  evolves  into  the  transmitted  ray,  both  of  which  lie  near  the 
critical  angle.  In  addition,  the  totally  reflected  ray  above  the  interface  becomes 
the  evanescent  ray  in  the  bottom.  Of  course,  the  incident  field  must  be  added  to 
the  reflected  field  to  obtain  continuity  with  the  transmitted  field  at  the  interface. 


A  comparison  between  classical  ray  theory  and  the  exact  result  [obtained 
by  numerical  integation  of  Eq.  (2-74)]  for  the  transmitted  field  along  a  horizontal 
line  is  shown  in  Fig.  2-13.  As  in  Fig.  2-4  for  the  reflected  field,  the  geometric 
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Figure  2-13  Comparison  between  classical  ray  theory  and  the  exact  result  for 

THE  NORMALIZED  TRANSMITTED  FIELD.  VALUES  OF  Z,  AND  Zr  IN  WAVELENGTHS  A  ARE  SHOWN 
FOR  EACH  PLOT.  THE  INTERFACE  IS  THE  SAME  AS  IN  FlG.  2-4. 

spreading  and  phase  terms  predicted  by  classical  ray  theory  have  been  removed 
from  the  total  field,  so  that  only  the  effect  of  the  transmission  process  remains. 
The  geometric  parameters  zs  and  zT  for  each  plot  are  indicated  on  the  figure. 
Note  that  the  ratio  z,/ zr  =  4  is  held  constant  so  that  the  classical  ray  result  is 
the  same  for  each  geometric  configuration.  Classical  ray  theory  is  less  accurate 
at  lower  frequencies  (source  and  receiver  closer  to  the  interface  in  terms  of  wave¬ 
lengths)  because  the  evanescent  field  is  more  significant.  The  interference  evident 
at  the  longer  ranges  (larger  0O)  is  between  the  transmitted  and  evanescent  rays. 
Attenuation  in  the  bottom  would  enhance  the  importance  of  the  evanescent  ray 
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in  comparison  with  the  transmitted  ray  because  the  latter  has  a  much  longer  path 
length  in  the  bottom. 

2-8.3  Geometric  interpretation  of  the  transmitted  field 

The  interpretation  in  terms  of  eigenrays  is  slightly  more  complicated  for 
the  transmitted  field  than  for  the  reflected  field.  In  Fig.  2-14  we  have  drawn  an 
eigenray  from  the  source  S  to  the  interface  at  A  at  angle  7  and  backward  from 
the  receiver  R  to  the  interface  at  B  at  the  angle  71  found  from  Snell’s  law  in 
Eq.  (2-71).  As  in  Fig.  2-8  for  the  reflection  case,  we  allow  for  the  possibility  of  a 
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Figure  2-14  Intercrktation  of  the  transmitted  field  in  terms  of  eigenrays. 

displacement  along  the  interface  by  not  assuming  that  points  A  and  B  coincide, 
and  we  define  as  A  the  resulting  ray  displacement  AB.  The  same  comments  made 
in  Secs.  2-7.1  and  2-7.3  regarding  the  possibility  of  complex  eigenray  angles  7 
apply  in  the  present  case. 
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To  find  an  expression  for  the  ray  displacement  A,  we  manipulate  the 
saddle  point  criterion  in  Eq.  (2-77).  Using  Snell’s  law  and  geometry,  we  can  write 
the  factor  (sin  9  +  q')  as 

/-in-.  sin^cos^  =  ~sin(7i  -  7)  =  -BEfBR  -BE 
\  n  cos  7j  j  cos  71  zr/BR  zr  ^ 


Substituting  for  sin(7  -  0O)  =  —CE/Ri  in  the  first  term  of  Eq.  (2-77),  we  obtain 


/'( 7)  = 


iTTE 

~rT 


iBE  (  1  IV'(7) 

R{  +  kRx  W{-i) 


(2-80) 


Multiplying  by  — z/21  and  recognizing  that  CE  —  BE  =  A  cos  7,  we  find 


A(7)  = 


1 

k  cos  7 


^(7) 

W(  7) 


(2-81) 


Thus,  the  ray  displacement  for  the  transmission  case  Eq.  (2-81)  is  obtained  from 
the  reflection  case  Eq.  (2-60)  by  simply  replacing  W  for  V. 

As  before,  the  eigenray  equation  is  obtained  by  substituting  A  into  the 
saddle  point  criterion.  Taking  our  cue  from  the  reflection  case,  we  multiply 
Eq.  (2-77)  by  iRi/ cos  7  and  expand  the  sin()  terms: 


iRi 

cos  7 


n  7) 


fllSin(7-fl0)  j  zr  sin(7l  —  7) 
cos  7  cos  7  cos  71 

/?i  tan  7  cos  90  —  R 1  sin  90  +  zT  tan  ji  —  zT  tan  7  -fi  A 


0  . 
(2-82) 


After  substituting  Rx  cos  90  =  z,  +  zr  and  i?,  sin  90  -  r  (see  Fig.  2-14),  we  obtain 
the  eigenray  equation 


r  =  z,  tan  7  +  A  +  zr  tan  71 

=  ri(7)  +  £(7)  +  r2(7)  .  (2-83) 

The  characteristics  of  the  transmitted  eigenray  are  obtained  from  the 
saddle  point  approximation  in  the  same  manner  as  for  the  reflected  eigenray: 


P(0o,  Rl,Zr)  = 


sin  7 


1 


iRifil) 


1U(  7)  exp[?'(4:3vA  +  A’Asin7  +  kiBR)} 


G< 


Gf, 


trans¬ 

mission 

loss 


ray  path  phase  <f> 


(2-84) 
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The  ray  path  phase  $  in  Eq.  (2-84)  is  obtained  as  follows: 

$(7)  =  kRi  cos(7  —  60)  -  kzr(cos  7  -  n  cos  71) 

=  k(^A  +  ~AC  +  CD)  -  kBRc os  7! (cos 7  -  nc os  71) 

=  k'tt+  &Asin7  +  &TT£[cos(7i  -  7)  -  cos  7,  cos  7]  +  k^TTR cos2  7, 

=  it^4  +  A:A  sin  7  +  kBR  sin  7  sin  71  +  ki  BR  cos2  7j 
=  kSA  +  fcA  sin  7  +  fcj  BR  .  (2-85) 

The  geometric  spreading  factors  in  Eq.  (2-84)  have  the  same  interpretation  as  for 
the  reflection  case  shown  in  Fig.  2-9.  In  terms  of  the  ray  path  lengths,  they  can 
be  written 


In  the  reflection  case,  the  ray  displacement  A  (at  a  lossless  interface)  turns 
out  to  be  real  for  totally  reflected  and  lateral  wave  eigenrays.  In  the  transmission 
case,  however,  the  displacement  A  turns  out  to  be  complex,  which  means  that 
the  simple  interpretation  of  a  real  displacement  upon  transmission,  as  depicted  in 
Fig.  2-14,  does  not  hold.  Using  the  case  illustrated  in  Fig.  2-11  as  an  example,  we 
find  that  A  is  largely  imaginary  for  transmitted  rays.  Its  magnitude  is  very  small 
at  short  ranges  (indicating  the  accuracy  of  classical  ray  theory  there)  but  larger  at 
longer  ranges.  For  example,  the  eigenray  for  the  receiver  at  0O  =  60°  and  zT  =  1A 
in  Fig.  2-11  has  angle  7  =  55°  —  i0.5°  and  displacement  A  =  (0.03  +  70. 28) A.  Recall 
that  classical  ray  theory  may  be  obtained  by  neglecting  the  displacement.  From 
Fig.  2-13  we  see  that  neglecting  the  displacement  for  the  receiver  at  60  =  60°  and 
zT  =  1A  results  in  an  error  of  0.6  dB.  Thus,  although  the  displacement  is  complex 
and  cannot  be  interpreted  physically,  we  must  include  it  in  the  mathematical 
treatment  of  the  transmitted  field.  Real  displacements  are  found  only  for  receivers 
at  long  ranges  on  the  interface  (zr  =  0),  where  A  is  the  same  as  for  the  reflected 
lateral  wave  ray  on  the  opposite  side  of  the  interface. 

For  the  saddle  points  in  Fig.  2-11  that  represent  the  evanescent  field 
(which  is  not  accounted  for  by  classical  ray  theory),  the  magnitude  of  A  is  on  the 
order  of  1.5-6A,  and  the  real  part  is  somewhat  larger  than  the  imaginary  part. 
By  Snell's  law  the  plane  waves  in  the  bottom  that  correspond  to  the  “evanescent 
ray"  saddle  points  are  inhomogeneous,  with  real  parts  close  to  90°.  The  fact  that 
the  eigenray  characteristics  have  significant  imaginary  parts  is  an  indication  that 
the  evanescent  field  in  the  bottom  cannot  be  represented  as  a  simple  ray  field. 
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2-9  CONCLUSIONS 

In  this  chapter  we  have  developed  a  systematic  method  for  finding  the 
reflected  and  transmitted  fields  due  to  a  point  source  in  the  presence  of  a  plane, 
penetrable  interface.  The  basic  technique  is  to  allow  the  plane  wave  reflection  and 
transmission  coefficients  to  influence  the  saddle  point  location  by  moving  them 
into  the  phase  function  /  before  taking  this  function’s  derivative  to  find  the  saddle 
points.  The  result  is  that  the  saddle  points  are,  in  general,  complex.  Complicated 
effects  of  the  reflection  and  transmission  processes,  such  as  the  reflected  lateral 
wave  field  in  the  water  and  the  evanescent  field  in  the  bottom,  are  now  included 
in  a  natural  way  in  the  saddle  point  formulation.  By  tracking  the  steepest  descent 
paths  in  the  complex  0-plane  and  comparing  them  with  the  required  contour  of 
integration,  one  can  determine  whether  multiple  eigenrays  exist  at  a  given  receiver. 
Numerical  integration  along  the  steepest  descent  paths  provides  a  robust  and 
frequency- independent  method  for  evaluating  the  eigenray  field  when  the  usual 
saddle  point  approximation  is  not  valid. 

An  attractive  feature  of  our  approach  is  that  each  saddle  point  can  be 
associated  with  an  eigenray  having  a  certain  path  of  propagation.  Characteristics 
such  as  ray  travel  time,  geometric  spreading  losses,  and  reflection  or  transmis¬ 
sion  loss  provide  one  with  physical  insight  into  the  problem.  In  certain  transition 
regions,  the  field  is  too  complicated  to  be  simply  expressed  as  a  single  ray  field 
or  even  the  superposition  of  several  ray  fields.  Examples  are  the  reflected  field 
in  the  region  of  the  critical  angle  and  the  transmitted  field  in  the  region  where 
the  transmitted  ray  is  first  accompanied  by  the  evanescent  ray.  In  such  regions 
steepest  descent  paths  are  curved  and  contributions  to  the  field  integral  come 
from  a  range  of  angles.  Since  the  saddle  point  approximation  is  not  valid,  nu¬ 
merical  integration  along  the  steepest  descent  path  must  be  used  to  evaluate  the 
field  strength.  Despite  the  fact  that  a  detailed  physical  interpretation  of  complex 
eigenray  characteristics  is  not  always  possible,  the  real  parts  of  the  eigenray  angles 
do  indicate  the  approximate  path  of  the  energy  represented  by  the  eigenray. 

The  overriding  motivation  for  this  work  has  been  to  develop  a  method 
for  solving  the  single  interface  problem  in  a  way  that  could  be  extended  to  more 
complex  environments.  In  the  following  chapter  we  will  apply  the  ray  method  to 
a  flat  waveguide  (the  Pekeris  model)  and  to  a  sloping  waveguide  (the  penetrable 
wedge  problem).  For  these  environments,  the  total  field  is  expressed  as  a  sum 
of  ray  fields,  each  of  which  can  be  evaluated  using  the  methods  outlined  in  this 
chapter. 


CHAPTER  3 


APPLICATION  OF  THE  RAY  METHOD  TO  FLAT  AND  SLOPING 

WAVEGUIDES 


3-1  INTRODUCTION 

In  the  previous  chapter,  a  method  was  developed  for  finding  the  field 
due  to  a  point  source  in  the  presence  of  a  single  plane,  penetrable  interface.  The 
method  involves  allowing  the  plane  wave  reflection  or  transmission  coefficient  to 
influence  the  location  of  the  (possibly  complex)  saddle  points,  as  well  as  their 
associated  steepest  descent  paths.  In  this  chapter  we  find  the  field  in  flat  and 
sloping  waveguides  by  using  the  same  saddle  point  method  to  evaluate  each  one 
of  the  ray  fields  that  contributes  to  the  total  field. 

The  so-called  Pekeris  waveguide  serves  as  a  simple  model  for  acoustic 
propagation  in  shallow  water  ocean  environments.  It  consists  of  a  flat,  homoge¬ 
neous  fluid  layer  (water)  that  is  bounded  below  by  a  homogeneous  fluid  halfspace 
of  higher  velocity  (sand,  for  example)  and  above  by  a  pressure-release  interface. 
Traditionally,  this  model  has  been  analyzed  using  normal  mode  theory,  beginning 
with  the  classic  work  of  Pekeris.17  Brekhovskikh’s  treatment  of  the  problem,  in 
Chapter  5  of  Ref.  5,  includes  both  the  ray  and  mode  approach.  Most  computa¬ 
tional  implementations  of  normal  mode  theory  include  only  the  discrete  modes, 
which  limits  the  validity  of  the  solution  to  long  ranges.  In  more  recent  years, 
closer  attention  has  been  paid  to  the  remainder  of  the  solution,  referred  to  as 
the  continuous  spectrum  or  the  branch  line  contribution.18  In  the  case  where  a 
mode  is  near  cutofF  (the  transition  from  the  discrete  spectrum  to  the  continuous 
spectrum  at  the  critical  angle),  normal  mode  theory  can  be  invalid  at  all  ranges. 
The  concepts  of  “virtual  modes”19  and  “pseudoresonances” 20  have  been  proposed 
as  means  for  correcting  the  theory  in  this  case. 

The  modeling  of  the  Pekeris  waveguide  in  terms  of  rays  with  beam  dis¬ 
placement  has  been  carried  out  by  Tindle,  et  al. 9,10,11  Although  this  ray  formu¬ 
lation  has  boon  successfully  applied  in  various  circumstances,  it  does  not  model 
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the  field  accurately  in  the  region  of  the  critical  angle.  In  this  region,  a  ray  theory 
caustic  is  formed,  and  correction  formulas  based  on  uniform  asymptotic  expan¬ 
sions  must  be  used.  However,  the  caustic  corrections  do  not  completely  solve  the 
problem  because  the  field  predicted  from  the  shadow-zone  formula  does  not  match 
up  well  with  the  partially  reflected  field  at  shorter  ranges.  As  can  be  seen  in  Fig.  1 
of  Ref.  12,  the  field  at  ranges  shorter  than  the  caustic  range  is  not  correctly  ac¬ 
counted  for  by  the  shadow-zone  formula  alone.  From  Fig.  1  of  Ref.  11,  we  see  that 
partially  reflected  rays  undergo  no  beam  displacement  and  are  abruptly  cut  off 
at  the  critical  angle.  Adding  the  partially  reflected  ray  field  to  the  shadow-zone 
field  only  corrects  the  field  at  short  ranges;  closer  to  the  critical  angle  the  field  is 
still  incorrect,  and  there  is  a  discontinuity  where  the  partially  reflected  rays  are 
cut  off  at  the  critical  angle.  The  reason  that  good  results  have  been  obtained 
using  the  usual  beam  displacement  formulation9,10,11  is  that  the  partially  reflected 
field,  corresponding  to  the  continuous  spectrum  in  mode  theory,  is  very  weak  at 
sufficiently  long  ranges  in  a  flat  waveguide. 

In  the  present  approach,  as  described  in  Chapter  2,  we  find  the  correct 
field  in  the  problematic  region  of  the  critical  angle  by  incorporating  the  entire 
reflection  coefficient  (rather  than  just  its  phase)  in  the  saddle  point  criterion, 
finding  complex  saddle  points,  and  numerically  integrating  along  the  paths  of 
steepest  descent.  The  field  calculated  in  this  way  matches  smoothly  with  the  field 
at  shorter  ranges  due  to  partially  reflected  rays,  as  well  as  with  the  field  at  longer 
ranges  due  to  totally  reflected  rays  plus  lateral  wave  rays.  For  lossless  interfaces, 
the  representation  of  the  field  at  long  ranges  is  equivalent  to  the  previous  beam 
displacement  formulation. 

The  second  environment  to  which  we  will  apply  our  ray  method  is  the 
sloping  waveguide  or  wedge,  which  is  a  useful  model  for  coastal  and  continental 
shelf  areas.  The  wedge  problem  has  received  considerable  attention  from  the 
underwater  acoustics  community  in  recent  years.  Jensen  and  Kuperman21  have 
used  a  parabolic  equation  model  to  demonstrate  the  cutoff  of  modes  in  upslope 
propagation.  Pierce22  and  Kamel  and  Felsen23  have  studied  the  same  issue  from  a 
theoretical  standpoint.  Arnold  and  Felsen24  discuss  the  field  in  a  wedge  in  terms 
of  rays  and  local  modes,  and  Arnold,  et  a/.,25,26,27  have  developed  the  theory  of 
intrinsic  modes  for  the  penetrable  wedge.  Xiang,  et  at ,,28  have  recently  reported 
progress  in  simplifying  the  intrinsic  mode  formulation. 

The  work  of  Tindle  and  Deane12  was  the  first  successful  application  of 
ray  theory  to  a  shallow  water  wedge.  Using  the  beam  displacement  formulation 
from  Ref.  10,  good  general  agreement  was  found  between  ray  theory  and  adiabatic 
normal  mode  theory  for  cases  where  the  adiabatic  approximation  is  valid  (small 
slopes  and  no  mode  cutoff).  In  such  cases,  the  problem  in  the  beam  displacement 
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formulation  in  the  area  of  the  critical  angle  apparently  did  not  affect  the  results 
significantly.  According  to  Sec.  1  of  Ref.  12,  partially  reflects!  energy  was  ac 
counted  for  by  way  of  the  shadow-zone  caustic  correction  fat  tors,  and  no  partially 
reflected  rays  were  included  in  the  calculations.  For  moderate  and  large  ocean 
bottom  slopes,  however,  proper  accounting  of  energy  incident  on  the  bottom  at 
angles  near  and  steeper  than  the  critical  angle  is  essential  because,  due  to  the 
change  in  ray  angle  upon  each  bottom  reflection,  many  rays  important  to  •In¬ 
field  calculation  are  reflected  from  the  bottom  at  such  angh-s.  In  contrast  to  the 
flat  waveguide,  these  partially  reflected  rays  are  important  at  all  ranges  from  tin- 
source  in  the  wedge,  not  just  at  short  ranges. 

The  first  three  sections  of  this  chapter  deal  wit  h  tin-  flat  waveguide.  Sec 
tion  3  2  describes  the  ray  method  for  calculating  the  field,  Section  3  3  contains 
comparisons  of  the  ray  model  with  the  SAFARI  model  and  a  normal  mode  model, 
and  Section  3—1  reviews  the  structure  of  the  eigen  rays  at  several  receivers  in  tin- 
waveguide.  The  next  three  sections  treat  the  sloping  waveguide.  Section  3  5  con 
tains  a  derivation  of  the  integral  expression  for  the  ray  fields,  Section  3  (i  contains 
a  comparison  of  the  ray  model  with  a  two-way  coupled  mode  model,  and  Sec 
tion  3-7  gives  examples  of  the  eigenray  structure.  Concluding  remarks  are  made 
in  Sec.  3-8. 

3-2  RAY  METHOD  FOR  THE  FLAT  WAVEGUIDE 

The  environment  assumed  is  the  standard  Pekeris  waveguide  illustrated 
in  Fig.  3-l(a).  The  parameters  needed  to  describe  it  are  the  ratio  of  sound  speeds 
n  =  c/ci,  the  ratio  of  densities  m  =  px/ p,  and  the  water  depth  h.  In  the  water, 
the  wave  number  k  is  related  to  the  wavelength  A,  the  frequency  /,  and  the  sound 
speed  c  by  k  =  2n  /  \  =  2irf/c.  The  source-receiver  geometry  is  specified  by  the 
horizontal  range  r  and  the  source  and  receiver  depths  z,  and  zr,  respectively.  All 
distances  are  expressed  in  terms  of  the  wavelength  A.  We  assume  that  n  <  1, 
which  is  typical  for  water- sand  interfaces  in  shallow  water  ocean  environments.  In 
such  a  case,  the  critical  angle  0cr  =  arcsin(rc)  is  real.  Attenuation  may  be  included 
in  either  medium  by  allowing  the  appropriate  sound  speed  to  have  an  imaginary 
part  (see  Sec.  2-5),  which  causes  n  to  be  complex  as  well. 

I  he  total  field  at  a  receiver  in  the  Pekeris  waveguide  is  expressed  as  a  sum 
of  ray  fields  due  to  the  source  and  its  images,  which  are  illustrated  in  Fig.  3-l(b). 
The  straight  lines  connecting  the  image  sources  with  the  receiver  represent  the  ray 
paths  of  classical  ray  theory.  We  specify  each  ray  path  by  the  number  of  bottom 
reflections  Nh,  the  number  of  surface  reflections  Ns,  and  the  ray’s  direction  of 
travel  at  the  receiver.  I  he  ray  field  due  to  each  image  source  is  represented  as 
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Figure  3-1  Geometry  for  the  flat  waveguide  (a)  Source-receiver  geometry 
(b)  System  of  image  sources,  each  of  which  represents  a  different  ray  path.  The 

THIN  DASHED  LINES  REPRESENT  IMAGED  SURFACES,  AND  THE  THIN  SOLID  LINES  REPRESENT 
IMAGED  OCEAN  BOTTOMS. 
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a  contour  integral  over  plane  waves  in  the  complex  0- plane,  where  0  is  the  plane 
wave  angle  of  incidence.  We  solve  the  integrals  in  terms  of  eigenrays ,  which  are 
specific  rays  at  specific  angles  that  travel  from  source  to  receiver.  It  is  important 
to  realize  that  there  may  be  more  than  one  eigenray  with  a  given  ray  path  (a 
given  sequence  of  reflections).  For  example,  we  saw  in  Chapter  2  that  a  lateral 
wave  eigenray  always  accompanies  a  totally  reflected  eigenray. 

The  expression  for  the  field  reflected  once  from  a  plane,  penetrable  in¬ 
terface  is  given  by  Eq.  (2-28).  To  modify  this  expression  for  a  ray  path  that  has 
Nt,  bottom  reflections,  with  plane  wave  reflection  coefficient  V,  and  N,  surface 
reflections,  with  reflection  coefficient  —1,  we  write 

P{Oo,Rl)  =  {-\)N-J—c'?[1-,™c'kR'r^e-0o)lV(0))N'V^OdO  ,  (3-1) 

V  27rr  J-f+.oo 

where  V  is  given  by  Eq.  (2-18).  For  each  of  the  image  sources  in  Fig.  3- 1(b)  the 
appropriate  values  of  the  incident  angle  0O  and  the  distance  R{  that  specify  the 
receiver  location  can  be  calculated  from  simple  geometry.  The  total  field  is  given 
by  an  infinite  sum  of  ray  fields  of  the  form  in  Eq.  (3-1)  due  to  the  infinite  number 
of  image  sources  that  exist.  In  practice,  the  sum  can  be  truncated  after  0O  has 
become  sufficiently  small  and  Nt,  sufficiently  large  that  the  ray  is  insignificant  due 
to  multiple  partial  reflections  at  the  bottom. 

As  in  Chapter  2  we  use  the  method  of  steepest  descent  to  solve  the  integral 
for  each  ray  field.  Before  applying  the  method  we  incorporate  the  reflection  process 
into  the  saddle  point  criterion  by  expressing 

[F(0)]"»  =  e">lnVW  ,  (3-2) 


so  that  we  may  write  the  ray  field  as 


p(0O,*l)  =  (-1  )W-Ji-C-f  r'°°  e*M.cos(*-M+^<nm]  v/~ 

V  2rrr  J-^+io o 

The  phase  function  in  the  steepest  descent  method  (see  Sec.  2-3)  is  now 

f(0)  =  lCos(0-0o)  +  -^r\nV(O)  , 


(3-3) 


kR, 


and  the  criterion  for  the  saddle  point  7  is 


m=-i*in(0-0o)  +  ^^=0  at  0  =  7  , 


(3-4) 


(3  5) 
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where  prime  indicates  differentiation  with  respect  to  9.  It  can  be  seen  by  compar¬ 
ing  Eq.  (3-5)  to  Eq.  (2-48)  that  it  is  no  more  difficult  to  find  the  saddle  points 
for  the  multiply  reflected  field  than  for  the  singly  reflected  field  since  the  only 
difference  is  the  factor  Nb- 

As  in  Sec.  2-7.1,  Eq.  (3-5)  can  be  interpreted  geometrically  as  an  eigenray 
equation,  where  a  horizontal  displacement 


A(7)  = 


i 

k  cos  7 


Y1 7) 

^(7) 


(3-6) 


occurs  at  each  of  the  Nb  bottom  reflections.  The  structure  of  the  saddle  points  as 
a  function  of  horizontal  receiver  range  r  is  similar  to  that  for  the  singly  reflected 
field  illustrated  in  Figs.  2-5  and  2-7. 

Once  the  saddle  point  7  is  found,  the  integral  in  Eq.  (3-3)  is  evaluated 
in  one  of  two  ways,  as  detailed  in  Chapter  2.  The  first  method  is  to  apply  the 
first-order  saddle  point  approximation  given  in  Eq.  (2-3S)  to  Eq.  (3-3).  The  result 
is 

=  (-l)"-  iv(-»))w‘  (.■*»■"»<•>-*>  ,  (3-7) 

where  f"  is  obtained  by  taking  an  additional  derivative  of  Eq.  (3-5).  The  second 
method,  which  is  required  when  the  saddle  point  is  near  the  critical  angle,  is  to  nu¬ 
merically  integrate  along  the  path  of  steepest  descent  in  the  complex  0-plane.  As 
discussed  in  Sec.  2-7.2,  the  \/sin  9  term  in  the  integrand  of  Eq.  (3-3)  is  evaluated 
at  the  saddle  point  7  and  taken  outside  the  integral  before  the  integration  is  per¬ 
formed.  We  do  so  because  the  integral  over  9  represents  the  geometrical  spreading 
Ge  of  the  ray  bundle  in  the  plane  of  Fig.  3-1,  while  the  factor  G 4,  =  \J sin  7 /r  rep¬ 
resents  the  spreading  in  the  perpendicular  plane.  We  will  see  that  the  separation 
of  the  spreading  factors  in  this  way  is  useful  when  constructing  the  integral  for 
the  sloping  waveguide  in  Sec.  3-5. 

The  algorithm  for  finding  the  field  at  a  series  of  receivers  at  constant 
depth  zT  proceeds  as  follows.  First,  the  field  due  to  the  direct  path  arrival  is 
found  at  all  the  receivers.  The  amplitude  of  this  field  serves  as  a  reference  for 
deciding  when  to  neglect  rays  that  are  significantly  weaker.  Next,  the  algorithm 
loops  through  the  image  sources  [see  Fig.  3-1  (b)],  finding  the  field  at  all  receivers 
for  each  ray  path.  The  ray  field  calculation  starts  at  the  receiver  for  which  0o  %  9C r, 
and  the  field  is  found  by  numerical  integration  along  the  steepest  descent  path. 
The  field  at  longer  ranges  (larger  angles  90  >  0cr)  is  found  next,  each  time  checking 
if  a  lateral  wave  ray  also  exists.  Once  the  lateral  wave  exists  also,  the  saddle  point 
approximation  is  compared  to  the  exact  (numerical)  result.  If  the  approximation  is 
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accurate  enough,  it  can  then  be  used  for  all  the  receivers  farther  out  in  range.  This 
procedure  for  testing  the  accuracy  of  the  saddle  point  approximation  is  applied 
to  the  lateral  wave  ray  as  well  as  to  the  totally  reflected  ray.  In  addition,  the 
strength  of  the  lateral  wave  ray  is  compared  to  the  reference  strength,  and  when 
the  lateral  wave  ray  is  sufficiently  weak,  it  can  be  neglected  for  receivers  at  longer 
ranges. 

The  same  type  of  procedure  is  used  for  the  partially  reflected  ray  field  at 
the  shorter  ranges,  for  which  0O  <  r-  As  A4  gets  larger,  the  partially  reflected  field 
drops  off  more  and  more  sharply  with  decreasing  range  r.  When  iVb  is  sufficiently 
large  that  the  ray  field  at  all  the  receivers  is  insignificant  due  to  numerous  partial 
reflections,  the  process  is  terminated  and  the  total  field  has  been  found. 

The  field  in  the  bottom  of  the  flat  waveguide  may  be  found  using  the 
same  approach.  Since  all  the  energy  arrives  from  above,  only  the  upper  half  of 
image  sources  in  Fig.  3-1  (b)  contributes.  From  the  expression  for  the  transmitted 
field  at  a  single  interface  Eq.  (2-75),  we  derive  the  following  for  the  flat  waveguide 
case: 

p(0o,Ri,zT)  =  (-l)iV*  \/r—  e'*  / 2  exp (fc/?,  i'cos(0-0o) 

V  ZtTT  J  —  y  +  ioo  L 

izr 

—  —(cos  0  —  n  cos  ) 

R\ 

+  j^(^lnVW  +  ]nW'W)  }de  .  (3-8) 

where  Nh  is  the  number  of  bottom  reflections  of  the  ray  path  (not  including 
the  final  transmission),  IF  is  the  plane  wave  transmission  coefficient  given  in 
Eq.  (2-73),  zT  is  now  the  positive  distance  below  the  bottom  interface,  and  6y  is 
the  plane  wave  angle  in  the  bottom  given  by  Snell’s  law,  sin  0X  =  sin  6 In.  In  solving 
Eq.  (3-8)  in  terms  of  saddle  points,  we  always  find  a  saddle  point  corresponding 
to  the  usual  transmitted  ray.  When  0O  is  large  there  may  also  be  a  saddle  point 
that  accounts  for  the  evanescent  leakage  of  energy  into  the  bottom,  as  discussed 
in  Sec.  2-8. 

3-3  COMPARISON  OF  MODELS  FOR  THE  FLAT  WAVEGUIDE 

In  order  to  test  the  validity  of  the  ray  model,  which  we  refer  to  as  the 
ACCURAY  model,  we  have  compared  its  predictions  of  propagation  loss  with 
those  of  other  models,  namely  the  SAFARI  fast  field  program29  developed  at 
the  SACLANT  Centre  in  fia  Spezia,  Italy,  and  a  normal  mode  model  called 
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NEMESIS30  developed  at  the  Applied  Research  Laboratories,  The  University  of 
Texas  at  Austin.  The  SAFARI  model,  whose  algorithm  solves  the  Bessel  trans¬ 
form  of  the  wave  equation  in  range- independent  environments  by  direct  numerical 
integration,  should  be  accurate  at  all  frequencies  and  all  but  very  short  ranges. 
The  normal  mode  model  includes  the  discrete  spectrum  in  its  calculations,  but 
does  not  account  for  the  branch  line  integral. 

The  first  case  for  comparison  of  the  models  was  chosen  so  that  a  mode 
would  be  just  on  the  verge  of  cutoff.  The  bottom  parameters  are  n  =  0.9375 


Figure  3-2  Comparison  of  the  ray  model  (ACCURAY),  the  SAFARI  model,  and 
a  discrete  normal  mode  model  for  a  flat  waveguide  of  depth  h  =  3.6A.  Three 

TRAPPED  MODES  EXIST,  BUT  MODE  3  IS  VERY  CLOSE  TO  CUTOFF  THE  EXTREMELY  CLOSE 
AGREEMENT  BETWEEN  THE  ACCURAY  AND  SAFARI  MODELS  MAKES  THEIR  PLOTS  DIFFI¬ 
CULT  TO  DISTINGUISH. 

[6ct  =  69.6°)  and  m  =  1.25,  and  the  source  and  receiver  depths  are  zs  —  0.72A  and 
zT  —  0.60A.  Figure  3-2  shows  a  comparison  of  the  three  models  for  a  water  depth 
of  h  =  3.60A,  which  is  slightly  deeper  than  the  cutoff  depth  for  the  third  mode 
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hc 3  =  3.59A.  The  figure  shows  that  the  ray  model  and  the  SAFARI  model  agree 
so  closely  that  the  plots  can  barely  be  distinguished.  The  normal  mode  model, 
however,  is  typically  off  by  4  dB  at  the  peaks  and  15  dB  at  the  valleys  of  the  plot. 
The  reason  for  its  poor  performance  in  this  case  is  that  the  pole  corresponding  to 
mode  3  is  very  close  to  the  branch  point  singularity  at  the  critical  angle.  These 
results  demonstrate  that  SAFARI  and  the  ray  model  correctly  account  for  energy 
in  the  neighborhood  of  the  critical  angle,  while  normal  mode  models  with  only 
the  discrete  spectrum  do  not. 

The  second  case  for  comparison  is  the  same  as  the  first  except  the  water 
depth  has  been  lowered  to  h  =  3.5A,  so  that  only  the  first  two  normal  modes  are 
trapped.  It  can  be  seen  from  Fig.  3-3  that  there  is  still  near  perfect  agreement 


Figure  3-3  As  in  Fig.  3-2,  except  the  water  depth  has  been  reduced  to  h  =  3.5A. 
Mode  3  is  now  cut  ofe,  but  still  contributes  significantly  at  the  shorter  ranges. 

between  the  ACCURAY  and  SAFARI  models.  The  mode  model’s  performance  is 
somewhat  better  than  for  the  previous  case,  but  it  is  still  quite  poor,  especially  at 
the  shorter  ranges.  Although  not  trapped  in  this  case,  Mode  3  does,  nevertheless, 
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influence  the  field  quite  strongly  out  to  a  range  of  at  least  200A,  which  corresponds 
to  57  water  depths. 

Since  most  realistic  ocean  bottoms  possess  some  attenuation,  we  consider 
in  Fig.  3-4  the  same  case  as  in  Fig.  3-2  except  that  an  attenuation  of  0.5  dB/A  has 
been  included  in  the  bottom.  The  continued  close  agreement  between  the  ACCU- 


Figure  3-4  As  in  Fig.  3-2,  except  the  bottom  has  an  attenuation  of  0.5  dB/A 
The  importance  of  mode  3  decreases  with  range  due  to  the  attenuation. 

RAY  and  SAFARI  models  indicates  that  attenuation  has  been  properly  included 
in  the  ray  model  approach.*  The  important  difference  from  the  lossless  case  is  that 
the  normal  mode  model’s  results  are  somewhat  more  accurate,  especially  beyond 
a  range  of  140A,  or  39  water  depths.  The  improvement  is  due  to  the  fact  that 
mode  3,  whose  proximity  to  the  critical  angle  causes  the  poor  performance  in  the 
lossless  case,  is  now  significantly  attenuated  with  increasing  range. 

‘The  ACCURAY  and  SAFARI  models  both  account  for  attenuation  by  way  of  complex  sound 
speeds  or  wave  numbers. 
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In  Ref.  31  Mannseth  shows  comparisons  of  propagation  loss  in  shallow 
water  between  the  formulation  of  ray  theory  with  beam  displacement  given  in 
Ref.  10,  the  discrete  normal  mode  solution,  and  the  exact  solution  (obtained  by 
adding  the  branch  line  integral  to  the  discrete  modes).  Similar  conclusions  are 
made  regarding  the  inadequacy  of  the  discrete  mode  solution  when  a  mode  is 
near  cutoff.  In  addition,  it  is  shown  that  the  solution  using  rays  with  beam 
displacement  is  also  quite  poor  in  the  region  of  mode  cutoff.  For  cases  where  no 
modes  are  near  cutoff  the  ray  solution  is  better,  but  not  as  good  as  the  discrete 
mode  solution.  Thus,  the  method  presented  here  for  obtaining  a  ray  solution  is 
an  improvement  over  the  formulation  used  in  Refs.  9-11.  The  main  difference  is 
that  the  present  method  accounts  correctly  for  the  transition  from  total  to  partial 
reflection. 

Our  final  comparison  of  models  appears  in  Fig.  3-5,  where  the  water 


Figure  3-5  Comparison  of  the  ray  model  (ACCURAY)  and  the  SAFARI  model  for 
A  FLAT  WAVEGUIDE  OF  DEPTH  h  =  0.5A.  No  TRAPPED  MODES  ARE  PRESENT. 

depth  has  been  reduced  to  h  =  0.5A  and  the  maximum  range  considered  is  100A, 
or  200  water  depths.  The  source  and  receiver  depths  are  z,  —  0.2A  and  zT  =  0.4A, 
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respectively.  Since  no  trapped  modes  are  present  (hcX  =  0.72A),  the  discrete 
normal  mode  model  cannot  be  applied.  Agreement  between  the  ACCURAY  and 
SAFARI  models  is  still  very  good.  The  field  in  this  case  decays  much  more  rapidly 
with  range  than  in  the  other  cases  considered.  The  ray  solution  at  the  longest 
range  requires  the  inclusion  of  ray  paths  with  up  to  40  bottom  reflections,  resulting 
in  a  total  of  more  than  160  eigenrays.  Despite  the  fact  that  each  individual  ray  is 
quite  strong,  the  destructive  interference  between  rays  results  in  a  weak  field.  The 
reason  that  the  ACCURAY  solution  is  not  smooth  at  long  ranges  is  that  small 
errors  in  the  individual  ray  fields  cause  imperfect  interference  between  rays. 

3-4  EIGENRAY  STRUCTURE  IN  THE  FLAT  WAVEGUIDE 

Figure  3-6  shows  the  saddle  points  and  their  corresponding  steepest  de¬ 
scent  paths  in  the  complex  0-plane  for  receivers  at  ranges  of  50A  and  200A  for  the 
environment  modeled  in  Fig.  3-2.  Each  of  the  saddle  points,  indicated  by  the  bold 
dots,  represents  the  contribution  of  one  eigenray  to  the  total  field  at  the  receiver. 
The  physical  path  that  each  eigenray  takes  in  the  waveguide  is  indicated  in  the 
figure  by  the  number  of  bottom  bounces  Nb.  The  steepest  descent  paths  drawn  in 
Fig.  3-6  are  the  actual  ones  that  would  be  used  to  find  the  field  from  numerical 
integration.  The  magnitude  of  the  integrand  of  Eq.  (3-3)  is  weaker  by  a  factor 
of  0.005  at  the  end  points  of  the  steepest  descent  paths  than  at  the  saddle  points. 

The  two  eigenrays  at  the  largest  angles,  for  which  Nb  =  0,  are  the  direct 
and  surface- reflected  rays.  Their  fields  are  simply  calculated  according  to  the 
spherical  wave  expression  elkRx  J Rx.  For  Nb  >  0  there  are  four  eigenrays  with  the 
same  number  of  bottom  interactions,  and  as  Nb  increases  the  eigenrays  become 
steeper  (their  angles  of  incidence  become  smaller).  For  each  totally  reflected  ray 
with  a  given  ray  path,  there  is  a  lateral  wave  ray  with  the  same  path  whose 
saddle  point  lies  next  to  the  critical  angle.  As  explained  in  Sec.  2-5,  its  steepest 
descent  path  loops  around  the  critical  angle  and  serves  to  complete  the  contour 
of  integration  for  the  integral  representing  that  ray  field.*  The  saddle  points  at 
smaller  angles,  whose  steepest  descent  paths  cross  the  branch  cut,  correspond  to 
the  partially  reflected  eigenrays.  Although  the  transition  from  total  to  partial 
reflection  is,  in  reality,  a  gradual  process,  we  use  the  criterion  of  steepest  descent 
paths  crossing,  or  not  crossing,  the  branch  cut  to  distinguish  between  partially 
and  totally  reflected  rays. 

The  eigenray  structure  for  the  receiver  at  the  short  range  of  50A  (14 

"To  avoid  confusion,  only  one  of  the  four  lateral  wave  paths  has  been  drawn  in  Fig.  3  6(a), 
and  only  three  of  the  28  lateral  wave  paths  have  been  drawn  in  Fig.  3  6(b). 
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8,  -d»g  (a)  r  *  50  X 


Figure  3-6  Saddle  points  and  steepest  descent  paths  for  the  eigenrays  at  re¬ 
ceiver  ranges  (a)  r  =  50A  and  (b)  r  =  200A  in  the  waveguide  modeled  in  Fig.  3-2 
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water  depths)  in  Fig.  3-6(a)  indicates  that  relatively  few  eigenrays  make  up  the 
total  field.  Only  those  with  Nb  =  1  are  totally  reflected.  Those  with  Nb  =  2  are 
partially  reflected,  but  their  strengths  are  only  1-6  dB  weaker  than  the  totally 
reflected  rays.  Additional  reflections  at  subcritical  angles  introduce  much  more 
loss:  eigenrays  with  Nb  =  2  and  Nb  =  3  are  approximately  24  dB  and  46  dB 
weaker,  respectively,  than  the  totally  reflected  rays.  The  lateral  wave  rays  (with 
./Vj=l)  are  20-23  dB  weaker  than  the  totally  reflected  rays.  Neglect  of  these  rays 
would  change  the  total  propagation  loss  by  about  0.8  dB. 

For  the  receiver  at  the  longer  range  of  200A  (56  water  depths)  in 
Fig.  3-6(b),  eigenrays  with  up  to  ten  bottom  interactions  are  shown.  Those  with 
Nb  >  8  are  partially  reflected  at  the  bottom.  With  no  attenuation  in  the  bottom, 
the  strengths  of  all  the  totally  reflected  eigenrays  are  essentially  the  same.  The 
total  horizontal  (beam)  displacement  at  the  bottom  of  these  totally  reflected  rays 
varies  from  24A  to  35A.  The  ray  displacement  per  reflection  decreases  from  24A  for 
Nb  —  1  to  5A  for  Nb  =  7.  The  lateral  wave  rays  with  Nb  =  1  are  39  dB  down  from 
the  totally  reflected  rays,  while  those  with  Nb  =  7  are  6-8  dB  down.  When  an 
attenuation  of  0.5  dB/A  is  included  in  the  bottom,  as  in  Fig.  3  4,  the  lateral  wave 
rays  at  both  ranges  are  severely  attenuated.  For  r  =  200A  the  totally  reflected 
rays  with  Nb  —  1  suffer  an  additional  attenuation  of  about  10  dB. 

On  both  of  the  plots  in  Fig.  3-6  the  angles  associated  with  the  trapped 
normal  modes  (see  Sec.  35.2  of  Ref.  5)  are  indicated  by  arrows.  A  comparison  of 
the  two  plots  points  out  the  basic  difference  between  the  ray  and  mode  descriptions 
of  the  field.  At  short  ranges,  few  rays  describe  the  field  accurately,  while  several 
discrete  modes  describe  the  field  only  approximately.  The  continuous  spectrum 
would  be  required  to  complete  the  mode  description.  At  long  ranges,  many  rays 
but  only  a  few  modes  are  required  to  describe  the  field.  The  mode  solution  would 
be  accurate  if  either  no  mode  were  near  cutoff  or  if  the  bottom  had  significant 
attenuation.  A  good  description  of  how  the  interference  between  rays  forms  a 
mode  may  be  found  in  Ref.  9. 


3-5  RAY  METHOD  FOR  THE  SLOPING  WAVEGUIDE 

3-5.1  Derivation  of  the  plane  wave  integrals 

We  now  consider  the  sloping,  or  wedge,  environment  illustrated  in 
Fig.  3-7(a).  Parameters  requited  to  describe  this  waveguide  are  n  =  c/ci, 
m  =  pi/p,  and  the  wedge  angle  a.  The  source- receiver  geometry  is  given  by 
the  vertical  water  depth  at  the  source  hs,  the  source  and  receiver  depths  z,  and 
zT ,  and  the  range  r.  The  image  sources,  which  lie  on  a  circle  centered  at  the  apex 
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Figure  3-7  Geometry  for  the  sloping  waveguide,  (a)  Source-receiver  geometry 
(b)  System  of  image  sources,  each  of  which  represents  a  different  ray  path.  The 

THIN  DASHED  LINES  REPRESENT  IMAGED  SURFACES,  AND  THE  THIN  SOLID  LINES  REPRESENT 
IMAGED  OCEAN  BOTTOMS 
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of  the  wedge,  are  shown  in  Fig.  3-7(b).  The  lines  connecting  the  image  sources 
to  the  receiver  again  represent  the  ray  paths  of  classical  ray  theory.  We  consider 
the  case  of  upslope  propagation  here,  but  the  downslope  case  poses  no  additional 
problems.  Diffraction  of  the  acoustic  energy  incident  directly  on  the  apex  of  the 
wedge  cannot  be  included  in  the  present  approach,  but  this  effect  is  negligible  in 
typical  ocean  environments,  where  o  is  small  and  the  bottom  is  penetrable. 

The  field  associated  with  each  image  source  (or,  equivalently,  ray  path) 
is  constructed  by  appropriately  modifying  the  plane  wave  integral  for  the  flat 
waveguide,  Eq.  (3-3).  For  a  given  image  source  in  Fig.  3  7(b),  the  distance 
and  the  angle  of  incidence  at  the  first  bottom  reflection  0O  can  be  computed  from 
simple  geometry.  The  main  difference  from  the  flat  waveguide  is  that  plane  waves 
undergo  successive  bottom  reflections  at  incident  angles  2o  steeper  (smaller)  than 
the  previous  one,  as  illustrated  by  the  ray  path  drawn  in  Fig.  3  7(b).  For  the 
general  ray  path  with  Ns  surface  reflections  and  A'/,  bottom  reflections,  we  may 
write  the  pressure  field  as 


p(90,Ri)  =  G*  J J„+  exp  j A-/? i  /  cos(0  -  0o) 

1  Nh  ,, 


(3-9) 


A  discussion  of  the  factor  G$  will  be  given  below.  Again,  there  may  be  more  than 
one  saddle  point  (eigenray)  when  we  solve  Eq.  (3-9). 

The  phase  function  of  the  integrand  in  Eq.  (3-9)  is 


f{9)  =  i  cos {9  -  0O )  +  7-5-  Y ln  V'[0  -  ~(J  ~  1  )Q1 

j-\ 


(3-10) 


and  the  saddle  point  criterion  is  now 


m  =  -<Sin(f - <w  +  jF E T$r| J -T^I  =  0  at,,  =  r  (3‘‘" 

Using  the  same  method  as  in  Sec.  2-7.1.  the  saddle  point  criterion  can  be  inter¬ 
preted  geometrically  as  an  eigenray  equation.  The  result  is  that  on  the  / 1 h  bottom 
reflection  from  the  source,  j  =  1,  ,V6,  the  eigenray  undergoes  a  displacement  par¬ 
allel  to  the  interface  of 


Mij)  = 


i 

kcas-tj  V(ij] 


where  ■)_,='>-  2 (j  - 


o 


(3  12) 
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Note  that  in  going  from  Eq.  (3-3)  to  Eq.  (3-9),  the  \/sin  6  factor  has 
been  taken  outside  the  integral  and  combined  with  the  \J  1/r  factor  to  get  (7^,  a 
geometric  spreading  factor  that  was  discussed  in  Sec.  3-2  for  the  flat  waveguide. 
In  Eq.  (2-70),  this  factor  was  derived  for  the  single  reflection  case  and  expressed 
as 

<7*  =  fi/L  ,  (3-13) 

where  L  was  related  to  the  ray  path  segments  between  source  and  receiver.  For 
the  sloping  waveguide,  L  is  more  complicated  due  to  the  change  in  direction  of 
the  ray.  It  takes  the  form 

Nb  M~i  ) 

L  =  ^(ray  path  lengths  in  water)  -F  — — —  .  (3-14) 

The  ray  field  in  the  bottom  of  the  wedge  can  be  derived  from  Eq.  (3-8) 
in  a  manner  similar  to  the  derivation  of  Eq.  (3  9).  The  result  is 

ik  *  f f 

p(0o,/?i,zr)  =  G+J  —  c'*  /  exps  kR\  icos{0-0o) 

Y  /7T  y_|+IOO  t 

iz  1  /Affc 

” -^(c°s(? -ncosi?1)  +  —  (J2ln  vl°  -  20  -  I)Q] 

+  In  W[0  -  2Nha\  )  J  dO  ,  (3-15) 

where  Nb  is  the  number  of  bottom  reflections  of  the  ray  path  (not  including 
the  final  transmission),  W  is  the  plane  wave  transmission  coefficient  given  in 
Eq.  (2-73),  zr  is  now  the  (positive)  perpendicular  distance  from  the  bottom  inter¬ 
face  to  the  receiver,  and  6\  is  the  plane  wave  angle  in  the  bottom  given  by  Snell’s 
law,  sin  =  sin  0/n.  The  geometric  spreading  factor  G $  for  this  case  again  takes 
the  form  given  in  Eq.  (3-13),  where  the  quantity  i  is  derived  from  Eq.  (2  87)  and 
is  given  by 


/Vf,+1  A(7  ) 

L  —  Y^(ray  path  lengths  in  water)  +  V'  - — 

%  s,n  7j 

ray  path  length  in  bottom 


(3  16) 


where  the  final  term  in  the  summation  over  j  is  to  be  evaluated  according  to 
Eq.  (2  81)  for  the  ray  displacement  upon  transmission. 
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3-5.2  The  mathematical  basis  for  the  multiplicity  of  lateral  waves 

An  interesting  aspect  of  propagation  in  a  sloping  environment  is  the  mul¬ 
tiplicity  of  lateral  waves,  caused  by  the  change  in  bottom  angle  upon  each  succes¬ 
sive  reflection.  Figure  3  8(a)  shows  the  three  possible  eigen  rays  for  the  ray  path 
with  N(,  =  2  and  Ns  =  1.  One  is  totally  reflected  at  both  bottom  reflections,  and 
the  other  two  are  lateral  waves,  each  of  which  strikes  the  bottom  near  the  critical 
angle  on  one  of  its  two  reflections.  Note  that  the  lateral  wave  at  the  first  bottom 
reflection  subsequently  undergoes  a  partial  reflection,  whereas  the  lateral  wave  on 
the  second  reflection  has  already  undergone  a  total  reflection.  Clearly,  tin*  latter 
is  stronger  under  normal  circumstances. 

Mathematically,  these  lateral  waves  are  accounted  for  by  keeping  track  of 
the  Riemann  sheets  on  which  l  is  calculated  for  the  various  hot  tom  reflect  ions.  As 
described  in  Sec.  2-4.2,  the  square  root  term  v»2  —  sin2  0  in  V  [see  Kq.  (2  IS)]  may 
be  evaluated  such  that  its  imaginary  part  is  either  positive  (the  upper  Riemann 
sheet)  or  negative  (the  lower  Riemann  sheet).  The  two  sheets  are  connected  along 
branch  cuts  that  begin  at  the  branch  points  at  ±0,.r.  as  shown  in  Fig.  2  3.  Kvery 
time  an  integration  path  crosses  a  branch  cut  in  the  complex  0-plane,  one  must 
change  Riemann  sheets  in  order  to  preserve  continuity  of  the  function  V(0). 

In  a  sloping  environment  the  reflection  coellicient  V'  is  evaluated  at  differ¬ 
ent  angles  for  each  bottom  reflection.  For  the  ray  path  with  A \  =  2,  l  is  evaluated 
at  0  on  the  first  reflection,  giving  rise  to  branch  points  at  ±0or,  and  at  0  —  2o 
on  the  second  reflection,  giving  rise  to  branch  points  at  ±0cr  -f  2n.  A  choice  of 
Riemann  sheet  must  be  made  for  each  of  these  two  evaluations.  Figure  3  8(b) 
shows  the  saddle  points  and  steepest  descent  paths  in  the  complex  0- plane  for  the 
eigenrays  traced  in  Fig.  3-S(a).  In  the  figure,  the  axes  representing  the  angle  at 
the  second  reflection  are  simply  shifted  to  the  right  by  an  amount  2o  compared 
to  the  axes  representing  the  angle  at  the  first  reflection.  They  have  also  b<vn 
artificially  shifted  upward  slightly  to  allow  one  to  distinguish  the  two  overlapping 
branch  cuts. 

The  Riemann  sheets  on  which  various  portions  of  the  steepest  descent 
paths  in  Fig.  3  8(b)  lie  are  indicated  by  the  plus  and  minus  signs  inside  paren¬ 
theses.  For  example.  (+  )  means  that  evaluation  of  \  is  on  the  upper  Riemann 

sheet  for  the  first  reflection.  V  (0).  and  on  the  lower  Riemann  sheet  for  l ho  second 
reflection,  \  (0  —  2 cv ) .  The  contour  of  integration  in  Kq.  (3  <))  requires  that  the 
path  begin  and  end  on  all  positive1  shevts  (f  +  ).  Keeping  in  mind  that  a  ehange- 
in  the’  appropriate  Rie’mann  shevt  e>e  e  urs  whe’uever  a  path  cretsse's  a  branch  eut .  it 
can  be  seen  that  by  aehling  the-  two  lateral  wave-  rays  1e>  the'  totaliv  re-lleeteei  rav. 
the  correct  contour  e>f  integration  is  comple'te'd.  Ne>te-  t  It;,  ill  three  saddle1  points 
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Figure  3  8  Multiple  lateral  waves  in  a  sloping  waveguide:  (a)  Ray  trace  of  hie 
TOTALLY  REFLECTED  RAY  PLUS  THE  TWO  LATERAL  WAVE  RAYS,  (b)  SADDLE  POINTS  AND 
STEEPEST  DESCENT  PATH S  IN  THE  COMPLEX  0-PLANE 
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lie  on  both  positive  sheets  (+  +).  In  the  general  case,  there  may  be  as  many  as 
Nt,  lateral  wave  rays  plus  one  totally  reflected  ray  for  a  ray  path  with  Nb  bottom 
reflections. 

The  algorithm  for  evaluating  the  total  field  at  a  series  of  receivers  along 
a  horizontal  line  in  a  wedge  is  similar  to  that  for  the  flat  waveguide  outlined  in 
Sec.  3-2.  In  the  present  case,  evaluation  of  the  field  for  each  ray  path  begins  at 
the  receiver  for  which  60  «  Qa  -f-  2 (Nb  -  1)0.  The  ray  at  this  angle  has  total 
reflections  at  all  interactions  except  the  last,  at  which  the  interaction  is  near  the 
critical  angle.  We  use  the  same  procedure  of  finding  the  field  at  longer  ranges 
first,  then  shorter  ranges.  At  each  successive  range,  we  check  for  the  validity  of 
the  saddle  point  approximation  and  for  the  weakness  of  the  eigenray.  Presence  of 
the  lateral  wave  rays  is  detected  by  checking  whether  the  steepest  descent  path 
crosses  the  appropriate  branch  cut. 

3-6  COMPARISON  OF  MODELS  FOR  THE  SLOPING  WAVE¬ 
GUIDE 

The  application  of  the  ray  method  to  acoustic  propagation  in  sloping 
waveguides  is  the  subject  of  Chapter  4,  in  which  we  compare  the  results  predicted 
by  the  ray  model  with  those  predicted  by  a  two-way  coupled  mode  model.  Here, 
we  show  just  one  of  the  model  comparisons.  The  problems  considered  are  the 
benchmark  problems  proposed  as  a  result  of  a  special  session  for  Underwater 
Acoustics  at  the  112th  Meeting  of  the  Acoustical  Society  of  America  in  Anaheim, 
California,  8-12  December,  1986.32  The  environment  considered  is  a  wedge  of 
angle  a  =  2.86°,  with  n  =  0.882,  m  =  1.5,  h,  =  3.33A,  and  z,  =  1.67A.  We  show 
in  Figure  3-9  a  comparison  for  a  case  of  upslope  propagation,  where  zT  =  0.5A  and 
r  =  0-67A.  The  agreement  is  almost  perfect,  even  for  receivers  in  the  bottom  at 
ranges  beyond  56.67  A.  The  close  agreement  is  stong  evidence  that  the  theoretical 
approach  outlined  here  is  valid  and  has  been  correctly  implemented. 

3-7  EIGENRAY  STRUCTURE  IN  THE  SLOPING  WAVEGUIDE 

Figure  3-10  illustrates  the  saddle  point  structure  for  a  case  of  upslope 
propagation  in  a  wedge  of  angle  a  =  2°,  with  n  =  0.866  (6CI  =  60°)  and  m  =  1.67. 
The  geometric  parameters  are  h,  =  5A,  z,  =  1A,  and  zT  =  0.5A.  Part  (a)  of  the 
figure  shows  the  saddle  points  and  their  steepest  descent  paths  in  the  complex 
0-plane  for  the  particular  ray  path  with  Nb  =  3  bottom  reflections  and  N,  =  2 
surface  reflections.  The  receiver  range  is  varied  from  r  =  30A  to  r  =  60A  in 
increments  of  2.5A,  as  labeled.  Recall  that  0  is  the  plane  wave  angle  at  the  first 
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Figure  3-9  Comparison  of  the  ray  model  (ACCURAY)  and  a  two-way  coupled 
MODE  MODEL  FOR  UPSLOPE  PROPAGATION  IN  A  WEDGE  OF  ANGLE  a  =  2.86°.  THE  RECEIVER 
AT  THE  MAXIMUM  RANGE  OF  67  X  LIES  IN  THE  BOTTOM  DIRECTLY  BELOW  THE  APEX  OF  THE 
WEDGE. 
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bottom  reflection,  so  the  critical  angle  lies  at  =  60°  for  the  first  reflection, 
at  da  +  2a  =  64°  for  the  second  reflection,  and  at  +  4a  =  68°  for  the  third 
reflection. 

From  Fig.  3-10(a)  we  can  deduce  that,  as  a  function  of  range,  the  reflected 
(non-lateral  wave)  eigenray  has  all  (three)  partial  reflections  for  r  <  40A,  two 
partial  reflections  and  one  total  reflection  beginning  at  r  =  42. 5A,  one  partial 
reflection  and  two  total  reflections  beginning  at  r  =  50A,  and  all  total  reflections 
for  r  >  57. 5A.  This  implies  that  there  are  lateral  waves  at  the  first  bottom 
reflection  for  ranges  r  >  42.5A,  at  the  second  bottom  reflection  for  r  >  50A,  and 
at  the  third  bottom  reflection  for  r  >  57.5A.  Since  the  steepest  descent  paths  of 
the  lateral  wave  rays  tend  to  overlap,  we  have  drawn  only  one  for  each  of  the  three 
bottom  reflections. 

Part  (b)  of  Fig.  3-10  shows  the  saddle  points  at  the  particular  range 
r  =  55A.  Saddle  points  that  correspond  to  eigenrays  with  zero  to  five  bottom 
interactions  are  indicated.  As  in  Fig.  3-6  for  the  flat  case,  there  are  groups  of 
fcur  saddle  points  for  each  value  of  Nb,  representing  the  four  multipath  eigenrays. 
As  Nt,  increases,  the  eigenray  incident  angles  become  smaller.  An  eigenray  with 
saddle  point  7  undergoes  its  jth  bottom  reflection  at  angle  7^  =  7  —  2 (j  —  1 )«,  so  it 
is  easy  to  determine  the  number  of  partial  reflections  a  given  eigenray  suffers.  For 
example,  all  four  eigenrays  with  Nb  =  2  experience  no  partial  reflections,  those 
with  Nb  =  3  experience  one,  and  those  with  Nb  =  4  experience  three.  Again,  only 
one  of  the  16  lateral  waves  at  0  =  60°  is  shown,  and  only  one  of  the  eight  lateral 
waves  at  0  —  64°  is  shown. 

This  phenomenon  of  eigenrays  with  larger  numbers  of  bottom  interact  ions 
gradually  ’  undergoing  more  partial  reflections  is  in  contrast  to  the  flat  waveguide 
case,  in  which  the  transition  from  rays  with  all  total  reflections  to  all  partial 
reflections  is  sudden.  An  important  consequence  is  that,  in  a  wedge,  eigenrays 
with  one  or  more  partial  reflections  are  significant,  even  at  long  ranges.  In  terms 
of  mode  theory  the  partially  reflected  contributions  correspond  to  modes  that  are 
beyond  cutoff  and  in  the  process  of  penetrating  into  the  bottom.  In  the  example 
of  Fig.  3-10(b),  the  rays  with  one  partial  reflection  (Nb  =  3)  are  just  2-4  dll 
weaker  than  the  totally  reflected  rays,  whereas  those  with  three  partial  reflections 
(Nb  =  4)  are  12-16  <111  weaker. 

3-8  CONCLUSIONS 

In  this  chapter  we  have  applied  the  complex  ray  methods  developed  in 
Chapter  2  to  the  flat  and  sloping  waveguides  common  in  underwater  acoustics 
applications.  For  the  flat  waveguide  (the  Fekeris  model)  agreement  between  the 
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ray  model  and  the  SAFARI  model  is  nearly  perfect.  Whereas  normal  mode  models 
that  calculate  only  the  discrete  spectrum  are  in  error  when  a  mode  is  near  cutoff, 
the  ray  model  is  valid  at  all  ranges.  The  ray  model  is  also  valid  when  the  water 
depth  is  as  low  as  half  a  wavelength. 

In  the  case  of  a  sloping  waveguide  we  have  derived  the  integral  expression 
for  each  ray  field,  and  have  shown  the  mathematical  origin  of  the  multiple  lateral 
waves.  In  contrast  to  the  flat  waveguide,  rays  in  a  wedge  change  angle  upon 
each  reflection.  Thus,  a  given  ray  that  is  totally  reflected  on  its  initial  bottom 
interactions  remains  significant  even  after  one  or  two  partial  reflections.  As  a 
result,  successful  modeling  of  a  wedge  environment,  even  at  long  ranges,  requires 
accurate  calculation  of  the  field  in  the  transition  region  from  total  to  partial 
reflection.  The  validity  of  the  ray  method  presented  here  has  been  demonstrated 
in  Sec.  3-6  by  its  close  agreement  with  an  independent  model. 

Although  the  environments  considered  in  this  chapter  are  simple  ideal¬ 
izations  of  the  real  ocean,  the  solutions  presented  here  in  terms  of  eigenrays  are 
valuable  because  they  provide  physical  pictures  of  the  propagation.  In  coim<v>t, 
models  that  numerically  solve  for  the  field  usually  provide  little  information  about 
the  actual  mechanisms  of  propagation.  Computationally  intensive  models,  such  as 
those  based  on  finite  element  methods,  the  parabolic  equation  approximation,  or 
coupled  mode  theory,  can  be  difficult  to  use  with  confidence  because  the  accuracy 
of  their  solutions  may  depend  on  the  initial  choice  of  mesh  size,  number  of  sample 
points,  number  of  integration  points,  or  number  of  modes  to  consider.  Thus,  a 
simple  model  like  the  ray  model  can  be  very  useful  in  verifying  the  accuracy  of 
these  more  versatile  and  powerful  models. 


CHAPTER  4 


RAY  MODEL  SOLUTIONS  TO  THE  BENCHMARK  WEDGE 

PROBLEMS 


4-1  INTRODUCTION 

In  this  chapter  we  present  in  more  detail  the  results  of  applying  the  rav 
method  developed  in  Chapters  2  and  3  to  the  benchmark  wedge  problems.  These 
results  were  first  reported  in  Ref.  4.  Benchmark  problems  were  proposed  at  a 
Special  Session  for  Underwater  Acoustics  at  the  112th  Meeting  of  the  Acoustical 
Society  of  America  in  Anaheim,  California,  8-12  December,  198G.32  Their  purpose 
was  to  assess  the  validity  and  accuracy  of  propagation  models  currently  in  use. 
In  the  case  of  the  penetrable-bottom  wedge,  an  analytical  solution  does  not  exist, 
so  only  by  comparing  approximate  solutions  obtained  from  independent  methods 
is  it  possible  to  verify  the  correct  solution.  In  Ref.  33  Jensen  and  Ferla  compare 
solutions  from  one-way  and  two-way  coupled  mode  models  and  from  two  models 
that  solve  the  parabolic  wave  equation.  These  models  are  quite  flexible  in  the 
environmental  geometries  to  which  they  can  be  applied,  but  they  provide  limited 
physical  insight  into  the  problem.  The  ray  approach  not  only  constitutes  an 
independent  method  for  solving  the  problem  but,  more  importantly,  provides  the 
user  with  a  physical  picture  of  the  propagation  paths  and  mechanisms. 

One  attractive  feature  of  the  ray  model  is  that  it  provides  detailed  in¬ 
formation  about  the  various  eigenrays  (rays  between  source  and  receiver)  that 
contribute  to  the  total  field.  The  relative  strengths  of  the  eigenrays  indicate 
which  paths  and  mechanisms  of  propagation  are  significant  and  which  arc  not. 
For  example,  from  the  eigenray  arrival  angles  one  can  determine  whether  the  phe¬ 
nomenon  of  backscatter  is  important.  In  a  wedge,  backscattered  eigenrays  have 
paths  that  travel  up  the  slope,  past  the  receiver,  and  then  back  down  the  slope 
before  striking  the  receiver.  This  propagation  in  the  "backward"  direction  is  rele¬ 
vant  because  some  models  are  intrinsically  unable  to  include  such  a  phenomenon. 
In  fact,  in  Ref.  33  the  difference  (of  2  dB  or  less)  between  the  one-wav  and  two- 
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way  versions  of  coupled  mode  theory  is  assumed  to  be  due  to  backscatter.  But  in 
Sec.  4-5  we  will  show  that  backscatter  is  definitely  not  a  factor  in  the  penetrable 
wedge  problems  proposed. 

Another  practical  aspect  of  wedge  propagation  addressed  by  the  ray 
model  is  the  importance  of  the  lateral  wave  field.  The  strengths  of  the  vari¬ 
ous  lateral  wave  contributions  to  the  field  at  a  particular  point  can  be  evaluated 
for  lossy  as  well  as  lossless  wedge  bottoms.  Besides  the  eigenray  strengths,  the 
ray  model  also  computes  eigenray  travel  times,  which  indicate  what  the  transient 
response  to  a  narrowband  pulse  would  be.  Since  the  lateral  wave  contributions 
arrive  earliest,  their  presence  may  be  noticeable  in  the  time  domain  even  if  their 
amplitudes  are  much  weaker  than  those  of  the  waterborne  eigenrays.  When  mea¬ 
suring  a  transient  signal  in  an  oceanic  wedge,  comparisons  of  lateral  wave  arrivals 
with  ray  model  simulations  might  be  used  to  determine  the  acoustic  parameters 
of  the  bottom. 

The  ray  model  differs  from  most  numerical  models  in  several  key  respects. 
First,  it  is  able  to  find  the  field  in  a  pointwise  manner;  t. e.,  it  is  not  based  on 
a  marching  algorithm  that  must  find  the  field  at  all  points  between  a  particular 
source  and  receiver.  Thus,  for  applications  where  the  field  at  only  a  handful  of 
points  is  desired,  the  ray  approach  is  likely  to  be  the  most  efficient.  Second,  the 
ray  model  is  more  efficient  at  higher  frequencies,  and  its  memory  requirements  are 
independent  of  frequency.  Most  numerical  models  require  more  computation  time 
at  higher  frequencies  and  are  limited  at  high  frequencies  by  the  storage  capacity 
of  the  particular  computer  on  which  they  are  implemented.  Further  details  on 
efficiency  considerations  will  be  given  in  Sec.  4-4. 

4-2  THE  BENCHMARK  WEDGE  PROBLEMS 

Three  wedge  problems  were  proposed  as  benchmarks.  The  first,  which  we 
will  not  discuss,  is  the  “ideal”  wedge,  which  has  a  pressure-release  bottom  as  well 
as  surface.  The  second  two  problems  represent  more  realistic  models  for  coastal 
and  continental  shelf  regions  of  the  ocean,  where  the  bottom  is  usually  sandy.  In 
the  second  benchmark  problem  the  bottom  is  penetrable  but  with  no  attenuation, 
while  in  the  third  problem  the  bottom  possesses  an  attenuation  of  0.5  dB/A. 
The  wedge  geometry  is  shown  in  Fig.  4-1.  The  wedge  angle  is  o  =  2.86°;  the 
frequency  is  25  Hz;  the  sound  speeds  in  the  water  and  bottom  are  c=1500  m/s 
and  Ci  =  1 700  m/s,  respectively,  such  that  n  =  0.88235  and  0„  =  61.93°;  and  the 
density  ratio  is  m  =  pxjp  —  1.5.  A  point  source  is  located  100  m  below  the 
surface  at  a  range  where  the  total  water  depth  is  200  m.  Propagation  loss  is  to  be 
calculated  in  the  upslope  direction  along  horizontal  lines  30  m  and  150  in  below 
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Figure  4-1  Geometry  for  the  benchmark  wedge  problems  proposed  by  the  Acous¬ 
tical  Society  of  America. 

the  water  surface.  The  horizontal  range  r  extends  from  zero  at  the  source  to  4  km 
directly  below  the  apex  of  the  wedge.  At  the  depth  of  30  m,  receivers  are  in  the 
water  for  r  <  3.4  km  and  in  the  bottom  for  r  >  3.4  km.  At  the  depth  of  150  m, 
receivers  enter  the  bottom  at  a  range  r=1.0  km. 

When  the  dimensions  of  the  wedge  are  expressed  in  terms  of  the  wave¬ 
length  A=60  m,  one  sees  that  this  is  indeed  a  low-frequency  environment.  The 
total  water  depth  decreases  from  just  3-1/3  A  at  the  source  to  zero  at  the  apex, 
and  the  total  horizontal  range  from  source  to  apex  is  66-2/3  A.  If  normal  modes 
are  calculated  at  each  range  as  though  the  water  depth  were  constant  (modes  for 
the  Pekeris  model),  then  at  the  source  there  would  be  just  three  trapped  modes. 
Mode  3  would  be  cut  off  at  a  range  of  r=0.81  km,  mode  2  would  be  cut  off  at 
r=2.09  km,  and  mode  1  would  be  cut  off  at  r=3.36  km. 

4-3  COMPARISONS  BETWEEN  MODELS 

At  the  time  the  benchmark  problems  were  proposed  (December  1986)  the 
following  propagation  models  that  could  be  applied  to  range-dependent  environ¬ 
ments  were  available:33  (1)  adiabatic  normal  mode  models,  (2)  coupled  normal 
mode  models  [one-way  and  two-way  versions],  (3)  models  based  on  the  parabolic- 
equation  (PE)  approximation  [narrow  or  wide-angle  formulations  and  split-step 
or  finite-difference  implementations],  (4)  finite-difference  solutions  to  the  full  wave 
equation,  and  (5)  finite-element  solutions  to  the  full  wave  equation.  Solutions  pro¬ 
vided  by  the  adiabatic  mode  models  are  not  expected  to  be  accurate  due  to  the 
simplifying  assumptions  on  which  the  models  are  based.  The  PE  models  should 
be  reasonably  accurate,  but  their  intrinsic  limitation  on  the  angular  spectrum  may 
introduce  error  in  a  sloping  environment.  The  finite-difference  and  finite-element 
methods  should,  theoretically,  provide  exact  solutions,  but  according  to  Ref.  33 
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the  existing  programs  were  too  computer  intensive  to  apply  to  these  problems. 
The  most  accurate  solution  available  is  from  a  two-way  coupled  mode  model, 
originally  developed  by  R.  B.  Evans,34,35  and  run  for  the  benchmark  problems  by 
Jensen  and  Ferla.33  There  is  also  a  one-way  version  of  the  coupled  mode  model 
which  is  more  efficient  but  which  ignores  propagation  in  the  backward  direction. 

In  this  section  we  compare  solutions  from  the  ray  model  to  solutions  from 
the  two-way  coupled  mode  model.  The  latter  set  of  solutions  was  kindly  provided 
to  the  author  by  Dr.  Finn  B.  Jensen  of  the  SACLANT  Undersea  Research  Center, 
La  Spezia,  Italy,  in  March  1988.  Comparisons  between  the  ray  model  and  the  two- 
way  coupled  mode  model  for  the  two  receiver  depths  appear  in  Fig.  4-2  for  the 
lossless  case  (Benchmark  2)  and  in  Fig.  4-3  for  the  lossy  case  (Benchmark  3). 
The  propagation  loss  in  Figs.  4-2  and  4-3  is  referenced  to  a  field  of  0  dB  at  1  m 
from  a  point  source  in  free  space.  The  overall  agreement  between  the  two  solutions 
is  excellent.  In  all  of  the  figures,  the  difference  is  less  than  1  dB.  The  agreement 
is  especially  good  at  receivers  in  the  water,  where  the  two  methods  differ  by  the 
line  thickness  or  less.  The  maximum  difference  of  1  dB  occurs  deep  in  the  bottom 
for  the  lossy  case  [Fig.  4— 3(b)].  Since  the  methods  used  by  these  two  models  are 
completely  different,  their  close  agreement  strongly  indicates  that  each  is  correct. 

To  calculate  the  field  at  the  400  points  in  each  of  the  figures,  the  ray  model 
requires  from  8  to  10  minutes  on  a  Cyber  180-830  computer.  Jensen  and  Ferla33 
report  that  the  two-way  coupled  mode  model  (run  on  the  FPS-164  array  processor) 
requires  24  hours  for  the  Benchmark  2  cases  and  8  hours  for  the  Benchmark  3  cases. 

4-4  ACCURACY  AND  EFFICIENCY  CONSIDERATIONS 

The  accuracy  of  the  ray  model  is  governed  by  two  separate  factors:  first, 
the  validity  of  the  approximations  that  were  introduced  into  the  theoretical  for¬ 
mulation  of  the  problem;  and  second,  the  accuracy  with  which  the  computer 
algorithm  calculates  the  ray  fields  given  by  the  approximate  integrals. 

As  discussed  in  Sec.  2-7.2,  two  approximations  were  introduced  into  the 
derivation  of  the  plane  wave  integral  Eq.  (3-9)  for  the  ray  field  in  a  wedge:  (1)  the 
substitution  of  the  asymptotic  form  of  the  Hankel  function  and  (2)  the  evaluation 
of  the  \/sin  0  term  at  the  saddle  point  7  and  its  subsequent  removal  from  inside 
the  integral  over  9.  The  ultimate  result  of  these  approximations  was  to  separate 
the  ray  spreading  terms  Go  and  G $  in  the  two  orthogonal  directions.  Since  prop¬ 
agation  straight  upslope  (or  downslope)  does  not  introduce  any  coupling  between 
the  ray  tube  spreadings  in  these  two  directions,  there  should  be  very  little  error 
in  this  formulation. 

As  mentioned  in  Sec.  3  5,  we  do  not  account  for  the  field  diffracted  from 
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Figure  4  2  Comparison  between  the  ray  model  and  the  two-way  cocpleh  mode 
model  eor  Benchmark  2.  (a)  30-vo  receiver  depth,  (b)  150-m  receiver  depth 
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Figure  4  3  Comparison  between  the  ray  model  and  the  two-way  coupled  mode 
MODEL  FOR  BENCHMARK  3:  (a)  30-m  RECEIVER  DEPTH,  (b)  150-111  RECEIVER  DEPTH. 
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the  apex  of  the  wedge.  This  field  is  produced  by  energy  that  travels  directly  from 
the  point  source  to  the  apex  and  is  diffracted  from  the  discontinuity  in  geometry 
and  boundary  conditions  there.  The  existence  of  the  diffracted  field  is  more  easily 
seen  when  considering  wedges  with  angles  greater  than  180°,  where  the  diffracted 
field  can  be  a  significant  part  of  the  total  field.  In  ocean  environments,  however, 
the  diffracted  field  should  be  insignificant  because  the  wedge  angles  of  interest  are 
typically  small  and  the  bottom  is  penetrable.  The  diffracted  field  should  not  be 
confused  with  the  backscatiered  field  that  consists  of  rays  that  travel  up  the  slope, 
“turn  around,”  and  travel  back  down  the  slope. 

The  accuracy  with  which  the  computer  algorithm  calculates  the  total 
field  can  be  controlled  by  the  user  of  the  ray  model.  For  each  receiver  the  pressure 
amplitude  of  the  direct  path  eigenray  is  used  as  a  reference,  and  any  eigenrays 
sufficiently  weaker  than  this  reference  are  discarded.  The  cutoff  amplitude  is  taken 
to  be  a  certain  number  of  dB  (specified  by  the  user)  below  the  reference.  For  the 
plots  in  Figs.  4-2  and  4-3,  this  “cutoff  margin”  was  set  to  a  very  conservative 
50  dB,  which  corresponds  to  ignoring  any  eigenray  whose  pressure  amplitude  is 
less  than  0.3%  of  the  reference.  Including  such  an  eigenray  would  change  the 
propagation  loss  by  just  0.027  dB.  In  practical  cases,  where  one  may  choose  to 
ignore  eigenrays  that  are  weaker  than  5%  of  the  reference  amplitude,  a  cutoff 
margin  of  26  dB  would  be  specified.  Including  such  an  eigenray  would  change  the 
propagation  loss  by  0.42  dB.  With  a  less  stringent  cutoff  margin,  fewer  eigenrays 
need  to  be  found,  and  the  computation  time  decreases  accordingly. 

The  time  required  to  find  the  field  at  a  single  receiver  depends  on  the 
number  of  significant  eigenrays,  which  in  turn  depends  on  the  receiver  location. 
At  short  ranges  in  the  wedge,  there  are  fewer  significant  eigenrays.  Those  with 
a  sufficiently  large  number  of  bottom  interactions  are  weak  because  they  strike 
the  bottom  at  steep  (small)  angles  where  only  partial  reflection  occurs.  At  longer 
ranges  these  eigenrays  are  significant  because  their  bottom  interactions  are  at 
shallower  (larger)  angles  where  total  reflection  occurs  with  little  or  no  loss.  For 
example,  the  field  at  a  range  of  0.5  km  in  Fig.  4-2(a)  requires  5  seconds  computa¬ 
tion  time  and  is  composed  of  contributions  from  just  14  rays,  having  up  to  three 
bottom  interactions.  Only  seven  of  these  rays  are  within  26  dB  of  the  strength  of 
the  direct  ray.  In  contrast  the  field  at  a  range  of  3  km  requires  almost  60  seconds 
computation  time  and  is  composed  of  contributions  from  71  rays  (having  up  to 
seven  bottom  reflections),  of  which  22  (having  up  to  five  bottom  reflections)  are 
within  26  dB  of  the  reference  strength. 

Another  factor  that  influences  the  computation  time  is  the  frequency.  At 
higher  frequencies  the  ray  algorithm  is  more  efficient  because  the  saddle  point 
approximation  to  the  integral  (rather  than  the  numerical  integration  along  the 
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path  of  steepest  descent)  can  be  used  more  often  to  evaluate  the  strengths  of  the 
eigenrays.  In  this  respect  the  ray  model  is  unique  since  models  based  on  normal 
mode  theory,  the  PE  approximation,  or  the  finite  element  method  can  become 
computationally  unfeasible  at  sufficiently  high  frequencies. 

Finding  the  field  at  multiple  receivers  along  a  horizontal  line  is  more 
efficient  than  finding  the  field  at  each  receiver  separately.  One  reason  is  that  an 
eigenray  found  for  one  receiver  can  serve  as  an  initial  guess  for  the  corresponding 
eigenray  at  a  nearby  receiver.  A  second  reason  is  that  in  order  to  use  the  more 
efficient  saddle  point  approximation  to  the  integral  for  the  eigenray  strength,  one 
must  first  verify  its  accuracy  by  comparing  it  with  the  result  from  numerical 
integration.  With  just  one  receiver  the  numerical  integration  method  must  be 
used  all  the  time.  A  final  reason  is  that  eigenrays  weaker  than  the  cutoff  level 
need  to  be  found  and  evaluated  before  they  can  be  recognized  as  being  too  weak. 
For  example,  in  the  field  calculation  at  the  single  range  of  3  km,  91  eigenrays  were 
actually  found,  even  though  only  71  were  within  50  dB  of  the  reference  strength. 
For  a  multiple-receiver  computation,  the  algorithm  recognizes  that  the  20  rays 
that  are  too  weak  do  not  need  to  be  found  at  the  other  receivers. 

Increased  efficiency  is  also  realized  by  interpolating  eigenrays  in  range. 
The  field  due  to  a  single  type  of  eigenray,  specified  by  its  number  of  surface 
and  bottom  reflections,  is  quite  smooth.  (It  is  only  when  the  fields  de*'  to  rays 
of  different  types  are  added  together  that  the  complicated  interference  patterns 
evident  in  Figs.  4-2  and  4-3  emerge.)  Consequently,  it  is  possible  to  find  eigenrays 
exactly  at  fairly  large  range  intervals  and  then  fill  in  the  intermediate  points 
by  interpolating  eigenray  characteristics  such  as  travel  time  and  geometric  loss. 
Interpolation  of  eigenrays  as  a  function  of  frequency  can  also  be  used  to  greatly 
reduce  the  computation  cost  of  simulating  the  propagation  of  broadband  pulses. 

Eigenray  interpolation  resulted  in  a  significant  improvement  in  efficiency 
when  f  igs.  4-2  and  4-3  were  produced.  Eigenrays  were  found  exactly  at  range 
intervals  of  120  m,  but  interpolation  was  used  to  narrow  the  interval  to  10  m. 
Thus,  the  “exact,”  field  was  computed  at  only  1/12  the  number  of  total  points. 
We  estimate  that  the  total  number  of  eigenrays  required  to  calculate  the  field  at 
the  400  points  in  Fig.  4  2(a)  is  on  the  order  of  20,000.  Using  the  interpolation 
scheme,  only  1700  of  these  need  to  be  found  and  evaluated. 


4-5  AN  \ LYSIS  OF  THE  FIELD  IN  THE  WEDGE 

An  attractive  feature  of  the  ray  model  is  that  it  provides  physical  insight 
into  how  acoustic  energy  travels  from  source  to  receiver.  For  each  receiver  position 
ill  the  wedge,  the  ray  model  compiles  a  list  of  eigenrays  that  contribute  to  the 
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total  field.  Mathematically,  each  eigenray  corresponds  to  a  saddle  point  in  the 
complex  0-plane.  Geometric  interpretation  of  the  quantities  associated  with  the 
saddle  point  yield  eigenray  characteristics  such  as  departure  and  arrival  angles, 
geometric  spreading  loss,  reflection  and  transmission  loss,  and  ray  travel  time  (see 
Sections  2-7  and  3-5).  In  this  section  several  issues  of  acoustic  propagation  in  a 
wedge  will  be  addressed  by  analyzing  the  field  in  terms  of  eigenrays. 

4-5.1  Characteristics  of  eigenrays  in  the  wedge 

As  an  example  of  the  eigenray  information  available,  the  eigenrays  for  the 
receiver  at  a  range  of  3  km  in  Fig.  4-2(a)  are  given  in  Tables  4-1  and  4  2.  Table  1  1 
lists  the  reflected  rays  while  Table  4-2  lists  the  lateral  wave  rays,  for  brevity,  only 
those  rays  with  six  or  fewer  bottom  reflections  are  included  in  Table  1  1.  Ravs 
with  more  reflections  are  at  least  26  dB  weaker  than  the  direct  ray.  In  Table  1  2, 
only  lateral  waves  with  four  or  fewer  bottom  reflections  are  included. 

One  interesting  observation  from  Table  4-1  is  that  rays  with  one  partial 
reflection  are  still  quite  strong  (rays  15-18  are  within  3-8  dB  of  the  direct  ray). 
This  indicates  that  partially  reflected  energy,  corresponding  to  untrapped  or  "vir¬ 
tual”  modes,  cannot  be  ignored  in  wedge  environments.  We  also  see  that  each 
additional  partial  reflection  results  in  a  loss  of  10-11  dB.  Another  observation  is 
that  the  eigenray  arrivals  are  quite  close  together  in  the  time  domain.  Rays  with 
up  to  six  bottom  reflections  arrive  within  four  periods  of  the  direct  arrival.  This 
is  an  indication  that  the  field  due  to  the  interference  between  many  ray  arrivals 
can  be  described  in  terms  rf  just  several  normal  modes. 

From  Table  4-2  we  see  that  there  are  a  very  large  number  of  lateral  wave 
ray  arrivals,  but  that  very  few  of  them  are  significant  compared  to  the  direct  rav. 
The  reason  for  the  large  number  is  that  a  lateral  wave  is  excited  whenever  energy 
is  incident  on  the  interface  at  angles  close  to  the  critical  angle  (0LT  =  62°  in  this 
case).  The  transmitted  energy  travels  in  the  bottom  parallel  to  the  interface  and 
continuously  reradiates  energy  back  into  the  water  at  an  angle  near  0cr.  In  a 
wedge,  a  great  number  of  lateral  waves  are  excited  by  energy  that  has  undergone 
different  numbers  of  surface  and  bottom  reflections  (see,  for  example,  Fig.  1  of 
Ref.  24).  Th  ere  may  be  as  many  as  A),  different  lateral  wave  eigenrays  having  .\\ 
bottom  interactions.  From  Table  4  2  we  see  that  the  strongest  lateral  wave  ray 
(ray  22)  is  about  14  dB  weaker  than  the  direct  ray.  ll  arrives  very  close  in  time 
to  the  reflected  ray  with  the  same  path  (ray  1  1  in  Table  4  '  ).  The  reflected  rav 
is  the  last  one  that  does  not  undergo  any  partial  reflections.  The  earliest  lateral 
wave  ray  (ray  1)  arrives  almost  five  periods  before  the  direct  arrival.  Although  it 
is  35  dB  weaker,  its  temporal  separation  from  the  stronger  arrivals  can  make  it 
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Table  4-1  List  of  the  reflected  eigenrays  and  their  characteristics  for  the  re¬ 
ceiver  AT  A  RANGE  OF  3  KM  IN  FlG.  4~2(a).  INITIAL  DIRECTION,  d, ,  SPECIFIES  WHETHER 
THE  RAY  LEAVES  THE  SOURCE  IN  AN  UPWARD  (U)  OR  DOWNWARD  (D)  DIRECTION.  ALSO 
GIVEN  ARE  THE  NUMBER  OF  SURFACE  AND  BOTTOM  REFLECTIONS,  Ar,  AND  N),\  THE  NUMBER 
OF  BOTTOM  INTERACTIONS  AT  WHICH  PARTIAL  REFLECTION  OCCURS,  Xp ;  THE  AMPLITUDE 
OF  THE  Ray  A  IN  dB;  THE  ANGLES  OF  INCIDENCE  AT  THE  FIRST  AND  LAST  BOTTOM  INTER¬ 
ACTIONS,  0j  AND  0I\  the  LOSSES  DUE  TO  GEOMETRIC  SPREADING,  Ag ,  AND  REFLECTION, 
Ar  \  AND  THE  TRAVEL  TIME  r,  EXPRESSED  IN  NUMBER  OF  PERIODS  OF  THE  25  Hz  WAVE  AND 
REFERENCED  TO  THE  ARRIVAL  OF  THE  DIRECT  EIGENRAY  . 


ray 

d, 

N. 

o, 

A 

4 

Ar 

T 

1 

U 

0 

0 

0 

- 

- 

-69.5 

-69.5 

0.0 

0.00 

2 

U 

1 

0 

0 

- 

- 

-69.6 

-69.6 

0.0 

0.03 

3 

D 

0 

1 

0 

86.6 

86.6 

-69.6 

-69.6 

0.0 

0.01 

4 

U 

1 

1 

0 

82.8 

82.8 

-69.6 

-69.6 

0.0 

0.09 

5 

D 

1 

1 

0 

85.4 

85.4 

-69.6 

-69.6 

0.0 

0.11 

6 

U 

2 

1 

0 

81.7 

81.7 

-69.7 

-69.7 

0.0 

0.26 

7 

D 

1 

2 

0 

83.7 

78.0 

-69.7 

-69.7 

0.0 

0.20 

8 

U 

2 

2 

0 

80.0 

74.3 

-69.8 

-69.8 

0.0 

0.39 

9 

D 

2 

2 

0 

82.6 

76.9 

-69.8 

-69.8 

0.0 

0.43 

10 

U 

3 

2 

0 

79.0 

73.3 

-69.8 

-69.8 

0.0 

0.68 

11 

D 

2 

3 

0 

81.0 

69.6 

-69.8 

-69.8 

0.0 

0.59 

12 

U 

3 

3 

0 

77.2 

65.8 

-69.6 

-69.6 

0.0 

0.88 

13 

D 

3 

3 

0 

80.0 

68.6 

-69.8 

-69.8 

0.0 

0.95 

14 

U 

4 

3 

0 

76.1 

64.7 

-69.6 

-69.6 

0.0 

1.29 

15 

D 

3 

4 

1 

78.5 

61.3 

-72.9 

-67.0 

-5.9 

1.20 

16 

U 

4 

4 

1 

75.3 

58.1 

-75.0 

-66.9 

-8.1 

1.61 

17 

D 

4 

4 

1 

77.8 

60.7 

-73.6 

-67.1 

-6.5 

1.69 

18 

U 

5 

4 

1 

74.3 

57.1 

-77.3 

-67.8 

-9.6 

2.14 

19 

D 

4 

5 

2 

77.3 

54.4 

81 .3 

-66.2 

-15.1 

2.05 

20 

U 

5 

5 

2 

73.9 

51.0 

-86.3 

-66.8 

-19.5 

2.54 

21 

D 

5 

5 

2 

76.5 

53.6 

-82.1 

-66.4 

-15.7 

2.63 

22 

V 

6 

5 

3 

73.3 

50.4 

-87. 4 

66.8 

20.6 

3.19 

23 

D 

5 

6 

3 

76.0 

47.4 

91 .7 

65.9 

-25.8 

3.06 

24 

li 

6 

6 

4 

72.9 

44.3 

97.7 

66.6 

31.2 

3.65 

25 

1) 

6 

6 

3 

75.2 

46.6 

93.7 

66.9 

26.8 

3.74 

26 

1' 

7 

6 

4 

72.4 

43.8 

98.8 

66.8 

32.0 

4.38 
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Table  4-2  List  of  the  lateral  wave  eigenrays  and  their  characteristics  for  the 
RECEIVER  AT  A  RANGE  OF  3  KM  IN  FlG.  4-2(a).  SEE  TABLE  4-1  FOR  AN  EXPLANATION  OF 
THE  E1GENRAY  CHARACTERISTICS. 


ray 

di 

N. 

Nk 

6i 

A 

A. 

Ar 

T 

1 

D 

0 

1 

0 

61.9 

61.9 

-104.0 

-104.0 

0.0 

-4.96 

2 

U 

1 

1 

0 

61.9 

61.9 

-102.0 

-102.0 

0.0 

-3.54 

3 

D 

1 

1 

0 

61.9 

61.9 

-103.5 

-103.5 

0.0 

-4.44 

4 

U 

2 

1 

0 

61.9 

61.9 

-101.4 

-101.4 

0.0 

-3.02 

5 

D 

1 

2 

0 

67.7 

61.9 

-98.9 

-98.9 

0.0 

-2.27 

6 

D 

1 

2 

1 

61.9 

56.2 

-109.6 

-103.2 

-6.4 

-4.07 

7 

U 

2 

2 

0 

67.7 

61.9 

-95.9 

-95.9 

0.0 

-1.16 

8 

U 

2 

2 

1 

61.9 

56.2 

-107.5 

-101.1 

-6.4 

-2.65 

9 

D 

2 

2 

0 

67.7 

61.9 

-98.2 

-98.2 

0.0 

-1.76 

10 

D 

2 

2 

1 

61.9 

56.2 

-109.1 

-102.7 

-6.4 

-3.47 

11 

U 

3 

2 

0 

67.7 

61.9 

-94.9 

-94.9 

0.0 

-.64 

12 

u 

3 

2 

1 

61.9 

56.2 

-106.9 

-100.5 

-6.4 

-2.05 

13 

D 

2 

3 

0 

73.4 

62.0 

-90.7 

-90.7 

0.0 

-.20 

14 

D 

2 

3 

1 

67.7 

56.2 

-104.2 

-97.8 

-6.4 

-1.38 

15 

D 

2 

3 

2 

61.9 

50.5 

-117.2 

-102.3 

-14.8 

-3.05 

16 

U 

3 

3 

0 

73.5 

62.1 

-85.2 

-85.2 

0.0 

0.60 

17 

U 

3 

3 

1 

67.7 

56.2 

-100.9 

-94.5 

-6.4 

-.27 

18 

u 

3 

3 

2 

61.9 

50.5 

-115.0 

-100.2 

-14.8 

-1.63 

19 

D 

3 

3 

0 

73.4 

62.0 

-89.4 

-89.4 

0.0 

0.31 

20 

D 

3 

3 

1 

67.7 

56.2 

-103.5 

-97.0 

-6.4 

-.79 

21 

D 

3 

3 

2 

61.9 

50.5 

116.7 

-101.9 

-14.8 

-2.37 

22 

U 

4 

3 

0 

73.6 

62.2 

83.5 

83.5 

0.0 

1.12 

23 

u 

4 

3 

1 

67.7 

56.2 

-99.9 

-93.5 

-6.4 

0.33 

24 

u 

4 

3 

2 

61.9 

50.5 

-114.4 

-99.6 

-14.8 

-.95 

25 

D 

3 

4 

1 

73.4 

56.3 

-95.3 

-89.0 

-6.3 

0.69 

26 

D 

3 

4 

2 

67.7 

50.5 

-111.4 

-96.6 

-14.8 

-.36 

27 

D 

3 

4 

3 

61.9 

44.8 

-125.9 

-101.6 

-24.3 

1.90 

28 

U 

4 

4 

1 

73.6 

56.5 

-89.5 

-83.8 

-5.7 

1.50 

29 

U 

4 

4 

2 

67.7 

50.5 

-107.8 

-93.0 

-14.8 

0.75 

30 

U 

4 

4 

3 

61.9 

44.8 

- 123.6 

99.2 

-24.3 

-.48 

31 

D 

4 

4 

1 

73.5 

56.3 

-94.0 

87.7 

-6.3 

1.29 

32 

D 

4 

4 

2 

67.7 

50.5 

-110.7 

95.9 

14.8 

0.32 

33 

D 

4 

4 

3 

61.9 

44.8 

-125.5 

101.1 

24.3 

1.15 

34 

U 

5 

4 

2 

67.7 

50.5 

106.9 

92.1 

14.8 

1.43 

35 

U 

5 

4 

3 

61.9 

44.8 

123.0 

98.7 

24.3 

0.27 

detectable  when  the  source  is  a  transient  pulse. 

Althouth  the  lateral  wave  eigenrays  are  weaker  in  the  presence  of  bottom 
attenuation,  it  is  still  necessary  to  include  them  to  obtain  an  accurate  solution. 
Along  the  interface,  lateral  wave  rays  travel  at  the  speed  of  sound  in  the  bottom 
and  are  attenuated  according  to  the  attenuation  in  the  bottom.  The  eigenray  list 
for  Benchmark  3  shows  that  ray  1  from  Table  4-2  suffers  an  additional  23  dB  loss 
due  to  the  attenuation  of  0.5  dB  per  wavelength  A  in  the  bottom.  This  translates 
to  a  path  length  along  the  interface  of  46A  or  2.76  km,  which  is  reasonable  for 
a  receiver  at  a  total  range  of  3  km.  The  strongest  lateral  wave  ray  in  Table  4  2 
(ray  22)  is  reduced  by  6  dB  as  a  result  of  attenuation,  which  makes  it  20  dB 
weaker  than  the  direct  ray.  Ignoring  such  a  ray  would  have  a  noticeable  effect  on 
the  total  propagation  loss. 

4-5.2  Eigenrays  in  the  bottom 

The  eigenrays  that  make  up  the  total  field  at  a  receiver  in  the  bottom 
are  listed  in  Table  4-3.  The  receiver  is  at  a  range  of  3.45  km  and  a  depth  of  30  m; 
it  lies  about  2.5  m,  or  0.042A,  below  the  interface.  Only  eigenrays  with  7  or  fewer 
bottom  reflections  are  included  in  the  table.  The  number  Nb  does  include  the  final 
bottom  interaction  at  which  transmission  (or  evanescent  leakage)  takes  place.  The 
four  types  of  rays  (see  Sec.  2-8.2)  listed  in  Table  4-3  are:  “TR”  for  transmitted, 
“LW”  for  lateral  wave,  and  “EV”  for  evanescent.  The  two  rays  labeled  ”EV2” 
represent  additional  evanescent  ray  saddle  points  and  are  described  in  Sec.  2-8.2. 

Some  general  characteristics  of  the  eigenrays  in  Table  4-3  will  now  be 
described.  In  upslope  propagation  the  evanescent  rays  have  total  reflections  at 
all  of  their  bottom  interactions.  At  a  fixed  receiver,  the  strongest  evanescent 
eigenray  (ray  20)  is  the  one  with  the  most  reflections  and  the  steepest  angle. 
Compared  to  the  transmitted  rays  with  the  same  ray  path,  the  evanescent  rays 
arrive  later  due  to  the  fact  that  more  of  their  paths  are  in  the  slower  medium. 
This  is  most  apparent  for  the  rays  that  travel  direct  from  the  source  into  the 
bottom:  evanescent  ray  2  arrives  6  periods  later  than  transmitted  ray  1.  For  the 
present  case,  where  the  receiver  is  very  close  to  the  interface,  the  evanescent  rays 
are  stronger  than  the  corresponding  transmitted  rays.  For  receivers  deeper  than 
several  wavelengths,  the  evanescent  rays  are  typically  negligible. 

The  lateral  wave  rays  listed  in  'Fable  4-3  arise  from  energy  that  undergoes 
a  lateral  wave  reflection  and  then  a  subsequent  transmission  into  the  bottom.  For 
example,  rays  9  and  12  are  reflected  as  lateral  waves  at  the  critical  angle  0cr  —  61.9° 
on  their  first  bottom  interaction  and  are  transmitted  into  the  bottom  on  their 
second  bottom  interaction  at  angle  flcr  -  2a  =  56.2°.  Rays  15  and  18  have  total 
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Table  4-3  List  of  eigenrays  and  their  characteristics  for  the  receiver  in  the 
BOTTOM  AT  A  RANGE  OF  3.45  KM  IN  FlG  4~2(a).  SEE  TABLE  4  1  FOR  AN  EXPLANATION 
OF  THE  EIGENRAY  CHARACTERISTICS. 


ray  typ. 


1 

TR 

D 

0 

1 

1 

61.9 

61.9 

-105.6 

-111.6 

6.0 

.00 

2 

EV 

D 

0 

1 

0 

86.6 

86.6 

-83.9 

77.5 

-6.4 

6.01 

3 

ev2 

D 

0 

1 

0 

91.2 

91.2 

94.7 

-82.0 

-12.7 

5.97 

4 

TR 

U 

i 

1 

1 

61.9 

61.9 

-104.1 

110.1 

6.0 

1.42 

5 

EV 

U 

i 

1 

0 

84.2 

84.2 

-76.8 

-73.6 

3.1 

6.04 

6 

ev2 

u 

i 

1 

0 

90.4 

90.4 

102.3 

-  86.1 

16.2 

5.80 

7 

TR 

D 

i 

2 

1 

67.6 

61.9 

-  101.7 

-107.7 

6.0 

2.69 

8 

EV 

D 

i 

2 

0 

85.9 

80.2 

-71.8 

73.0 

1.2 

6.07 

9 

LW 

D 

i 

2 

1 

61.9 

56.2 

-102.9 

106.3 

3.4 

.49 

10 

TR 

U 

2 

2 

1 

67.6 

61.9 

99.4 

-105.5 

6.0 

3.80 

11 

EV 

U 

2 

2 

0 

82.7 

77.0 

-69.7 

72.5 

2.8 

6.15 

12 

LW 

U 

2 

2 

1 

61.9 

56.2 

-  101.2 

-104.6 

3.4 

1.91 

13 

TR 

D 

2 

3 

1 

73.4 

61.9 

-95.7 

-101.7 

6.0 

1.75 

14 

EV 

D 

2 

3 

0 

84.2 

72.8 

-68.0 

-72.5 

4.6 

6.22 

15 

LW 

D 

2 

3 

1 

67.7 

56.2 

-98.6 

-101.9 

3.4 

3.17 

16 

TR 

U 

3 

3 

1 

73.4 

61.9 

-91.9 

-97.9 

6.0 

5.55 

17 

EV 

U 

3 

3 

0 

81.0 

69.6 

-66.9 

-72.2 

5.3 

6.36 

18 

LW 

U 

3 

3 

1 

67.7 

56.2 

-96.0 

-99.4 

3.4 

4.29 

19 

TR 

D 

3 

4 

1 

79.1 

62.0 

-85.2 

-91.2 

5.9 

6.18 

20 

EV 

D 

3 

4 

0 

82.5 

65.4 

-65.8 

-72.1 

6.3 

6.47 

21 

LW 

D 

3 

4 

1 

73.4 

56.2 

-92.0 

-95.4 

3.4 

5.24 

22 

TR 

U 

4 

4 

1 

79.0 

61.8 

-66.3 

-70.1 

3.9 

6.64 

23 

LW 

u 

4 

4 

1 

73.4 

56.3 

-87.7 

-91.2 

3.4 

6.04 

24 

TR 

D 

4 

5 

1 

81.0 

58.1 

-67.0 

-69.6 

2.6 

6.82 

25 

LW 

1 

4 

5 

1 

79.4 

56.5 

80.7 

-84.3 

3.7 

6.69 

26 

TR 

U 

5 

5 

2 

78.3 

55.4 

-71.0 

-68.1 

2.8 

7.10 

27 

LW 

u 

5 

5 

2 

73.4 

50.5 

94.0 

90.4 

3.5 

6,59 

28 

TR 

D 

5 

(i 

2 

80.1 

51.5 

75.3 

-68.2 

7.1 

7,31 

29 

LW 

D 

5 

6 

3 

73.4 

44.8 

106.3 

94.0 

12.3 

6.12 

30 

TR 

U 

6 

6 

3 

77.7 

49.1 

79.1 

-67.5 

1 1.6 

7.64 

31 

LW 

u 

6 

6 

3 

73.4 

44.8 

101.8 

89.7 

12.2 

7.22 

32 

TR 

n 

6 

7 

3 

79.6 

45.3 

84.7 

67.0 

17.8 

7.90 

33 

LW 

i) 

6 

7 

4 

73.4 

39.1 

1  15,1 

93.4 

21.9 

7.1  1 

34 

TR 

u 

7 

7 

1 

77.2 

42.9 

88.9 

67.2 

21.7 

8.26 

35 

LW 

i; 

7 

7 

1 

73.4 

39.1 

1  10.8 

89.0 

21.8 

7.91 
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reflections  at  6a  +  2a  =  67.7°  on  their  first  interaction  and  then  follow  the  type  of 
path  just  described.  Rays  21  and  23  undergo  two  total  reflections,  followed  by  the 
lateral  wave  reflection  and  subsequent  transmission,  while  ray  25  undergoes  three 
total  reflections  before  the  lateral  wave  reflection  and  subsequent  transmission.  As 
a  last  example,  ray  27  undergoes  two  total  reflections,  a  lateral  wave  reflection,  a 
partial  reflection,  and  then  a  transmission  into  the  bottom.  Most  of  the  lateral 
wave  rays  are  weak,  but  the  strongest  one  from  Table  4-3  (ray  25)  has  a  strength 
of  —80.7  dB,  which  is  just  15  dB  weaker  than  the  strongest  ray  on  the  list  (ray  20). 
The  multiplicity  of  lateral  waves  is  an  indication  of  the  complexity  of  the  field  in 
the  bottom. 

The  strongest  transmitted  (TR)  rays  listed  in  Table  4-3  are  those  having 
four  and  five  bottom  interactions.  Ray  22  is  within  0.5  dB  of  being  the  strongest 
ray  at  the  receiver.  Transmitted  rays  with  few  bottom  interactions  (rays  1  and  7, 
for  example)  strike  the  bottom  on  their  last  interaction  at  angles  near  0CT,  pen¬ 
etrate  at  very  shallow  angles,  and  suffer  much  loss  due  to  geometric  spreading 
( Ag ).  Those  with  many  bottom  interactions  (rays  30  and  34,  for  example)  suffer 
much  loss  due  to  multiple  partial  reflections  (Ar). 

4-5.3  Backscatter 

An  important  issue  in  the  penetrable  wedge  problem  is  the  phenomenon 
of  backscatter.  Backscattered  energy  corresponds  to  rays  that  leave  the  source 
traveling  upslope  and  eventually  strike  the  receiver  traveling  downslope.  The 
change  in  direction  occurs  as  the  result  of  multiple  reflections  off  the  sloping 
bottom.  A  backscattered  eigenray  is  easily  recognized  as  one  that  has  a  negative 
incident  angle  at  its  last  bottom  reflection  (0/  in  Tables  4-1  and  4-2).  A  brief 
analysis  of  the  eigenrays  at  various  ranges  indicates  that  backscatter  is  totally 
insignificant  for  the  penetrable  wedge  problems.  For  the  receiver  at  a  range  of 
3  km  in  Fig.  4-2(a),  the  steepest  eigenray  within  50  dB  of  the  reference  strength 
arrives  at  an  incident  angle  of  33°.  From  Table  4-1  we  see  that  the  steepest 
eigenray  within  29  dB  arrives  at  43.8°  (ray  26).  It  is  only  at  very  short  ranges 
(less  than  0.1  km)  that  eigenrays  traveling  vertically  are  significant. 

Using  the  image  method  it  is  easy  to  find  the  approximate  path  of  the 
eigenray  that  would  arrive  at  the  receiver  traveling  vertically.  In  Fig.  4-4  we  see 
that  the  image  sources  in  a  wedge  lie  on  a  circle  centered  at  the  apex.  The  ray 
arriving  vertically  at  the  receiver  is  imaged  backwards  to  the  closest  image  source 
on  the  circle.  From  geometry  this  image  source  lies  approximately  at  the  angle  / i 
given  by 


/f  =  arccos  (ItA/SA) 


(4-1) 
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The  approximate  number  of  bottom  reflections  experienced  by  the  first  backscat- 
tered  ray  is  then  Nt  =  /?/2a.  For  the  receiver  at  3  km,  /?  =  75°  and  Nt  =  13 
bottom  reflections.  Since  the  final  reflection  is  at  an  incident  angle  close  to  0° 
and  the  critical  angle  is  at  62°,  the  last  ten  reflections  are  at  angles  steeper  than 
critical  (62°  -f2a  «  10).  With  about  11  dB  loss  at  each  partial  reflection,  it  is 
clear  that  the  backscattered  eigenrays  are  negligible  for  the  present  case. 

Although  backscatter  is  not  a  factor  for  the  typical  water-sand  coastal 
wedge  specified  in  the  benchmark  problem,  the  previous  analysis  indicates  the 
conditions  under  which  backscatter  may  be  a  factor.  First,  in  a  wedge  of  larger 
angle  a,  fewer  reflections  are  required  to  make  a  ray  “turn  around.”  Second, 
a  wedge  with  a  higher  velocity  bottom  has  a  smaller  critical  angle  of  incidence, 
resulting  in  fewer  partial  reflections  at  the  bottom  interface.  And  finally,  backscat¬ 
tered  rays  experience  less  loss  when  the  wedge  bottom  has  less  attenuation. 

4-6  CONCLUSIONS 

The  excellent  agreement  between  the  ray  model  and  the  two-way  coupled 
mode  model  is  strong  evidence  that  both  are  capable  of  predicting  the  acoustic 
field  in  a  realistic  penetrable  wedge  with  an  error  that  is  typically  much  less  than 
1  dB.  The  ray  model  is  unique  in  that  it  expresses  the  total  field  at  a  point  as  a  sum 
of  ray  fields.  From  a  list  of  eigenrays  and  their  characteristics,  the  user  can  gain 
much  insight  into  the  relative  importance  of  the  various  paths  and  mechanisms  of 
propagation.  Unlike  most  other  models,  the  ray  model  is  more  efficient  at  higher 
frequencies  and  is  able  to  find  the  field  at  a  single  point  without  the  need  to  march 
through  the  entire  region  between  the  source  and  receiver. 

For  the  realistic  wedge  problems  that  were  proposed  as  benchmarks,  an 
eigenray  analysis  demonstrates  conclusively  that  backscatter  is  not  a  significant 
mechanism  of  propagation.  This  raises  a  basic  question  regarding  the  one-way  and 
two-way  versions  of  coupled  mode  theory.  The  intuitive  assumption  has  been  that 
the  one-way  version  neglects  energy  propagating  in  the  “backward”  direction,  so 
that  the  difference  between  the  two  versions  is  due  to  the  phenomenon  of  back¬ 
scatter.  However,  in  the  benchmark  wedge  problems  the  absence  of  backscatter 
calls  into  question  the  reason  for  the  difference  of  up  to  2  dB  between  the  one-way 
and  two-way  solutions  observed  in  Figs.  5  and  9  of  Ref.  33. 

The  two-way  coupled  mode  model  has  already  been  shown  to  be  in  good 
agreement  with  the  anaytical  solution  for  the  ideal  wedge  (one  having  perfectly 
reflecting  boundaries).36  However,  propagation  in  a  penetrable  wedge  is  quite  dif¬ 
ferent  in  that  the  reflection  process  at  the  bottom  is  very  complex.  For  example, 
the  entirely  new  mechanism  of  lateral  wave  propagation  is  introduced,  and  there 


exists  a  complicated  field  in  the  bottom  of  the  wedge  that  is  generated  not  only 
by  transmission  but  also  by  evanescent  leakage  of  totally  reflected  energy.  With 
such  a  profound  difference  between  the  two  problems,  it  is  not  clear  that  the  va¬ 
lidity  of  a  model  for  one  necessarily  guarantees  its  validity  for  the  other.  Thus, 
the  corroboration  of  the  two-way  coupled  mode  model  by  the  ray  model  for  the 
penetrable  wedge  is  an  important  result. 


CHAPTER  5 


THE  THREE-DIMENSIONAL  WEDGE  PROBLEM 


5-1  INTRODUCTION 

The  three-dimensional  wedge  problem  is  the  same  as  that  considered  in 
Sec.  3-5  and  Chapter  4  except  that  the  receiver  is  not  required  to  lie  straight 
upslope  or  downslope  from  the  so'^ce.  Detailed  theoretical  work  on  this  type 
of  geometry  has  dealt  mainly  with  the  special  case  of  perfectly  reflecting  surface 
and  bottom  boundary  conditions,  a  case  we  shall  call  the  “ideal  wedge.”  In  the 
works  of  Bradley37  and  Bradley  and  Hudimac,38  the  authors  provide  a  general 
review  of  the  literature  on  the  ideal  wedge,  going  back  to  Sommerfeld’s  classic- 
work.39  They  go  on  to  present  both  the  ray  and  mode  formulations  of  the  problem 
as  applied  to  underwater  acoustics  and  to  provide  approximate  mode  solutions 
for  various  source-receiver  configurations.  Buckingham40  derives  a  similar  normal 
mode  formulation  but  obtains  solutions  with  a  more  extensive  range  of  validity. 

In  Ref.  41,  Weston  studies  the  phenomenon  of  horizontal  refraction  in  a 
wedge,  where  ray  paths,  as  viewed  from  above,  are  bent  due  to  repeated  reflections 
at  the  wedge  surface  and  bottom.  He  derives  some  geometrical  characteristics  of 
ray  paths  in  a  wedge  and  qualitatively  analyzes  the  effects  of  a  more  realistic 
bottom  possessing  attenuation.  Harrison42  analyzes  shadow  zones  in  three  dimen¬ 
sional  environments  by  restricting  ray  elevation  angles  to  the  discrete  values  that 
correspond  to  the  angles  associated  with  the  normal  modes  at  the  source.  In  the 
work  of  Doolittle,  et  a/.,43  experimental  evidence  of  horizontal  refraction  is  given. 
Ray  theory  and  Buckingham’s  normal  mode  technique  are  used  to  analyze  the 
data  theoretically.  Although  the  existence  of  a  three-dimensional  model  based  on 
the  parabolic  equation  approximation44  was  noted,  the  model  was  not  available 
for  use. 

lu  this  chapter  we  apply  the  ray  methods  developed  earlier  to  find  the 
field  in  the  three-dimensional  wedge.  We  find  that  the  dependence  of  a  ray’s 
incident  angle  at  the  bottom  on  both  angular  variables  0  and  <p  prevents  the 
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integral  over  <j>  from  being  performed  analytically.  An  approximate  method  of 
solution  is  proposed. 

5-2  DERIVATION  OF  THE  PLANE  WAVE  INTEGRAL 

The  geometry  for  the  three-dimensional  wedge  problem  is  illustrated  in 
Fig.  5-1  (a).  The  z  =  0  plane  contains  the  ocean-air  interface,  the  y-axis  is  the 
shoreline  or  wedge  apex,  and  the  ocean  bottom  makes  an  angle  a  with  the  sur¬ 
face.  The  source  lies  in  the  y  =  0  plane  at  (x,,0,z,),  and  the  receiver  lies  at 
(xr,yr,zr).  The  system  of  imaged  sources  is  illustrated  in  Fig.  5—  1(b).  As  in 
the  two-dimensional  problem  of  Sec.  3-5,  the  images  lie  on  a  circle  centered  on 
the  wedge  apex,  and  each  line,  of  length  ft,,  connecting  an  image  source  to  the 
receiver  represents  a  ray  path  of  classical  ray  theory.  The  direction  of  the  line  is 
specified  by  the  angles  90  and  <£0,  obtained  by  translating  the  line  to  the  origin 
and  defining  the  usual  polar  angles  0  and  <f>  as  in  Fig.  5  1(a).  For  an  image*  source 
at  coordinates  (x,»,0,  z,»),  the  three  quantities  associated  with  the  ray  path  are 
computed  from  simple  geometry: 

Rl  -  \J(*r  -  X.')2  +  Vr  +  (^r  -  Za>)2  (5  la) 

Z  —  Z  I 

COS0O  -  ----  0  <  0O  <  (5-lb) 

til 

xr  —  xaf  yr 

cos  <f>o  =  -  sin  4>o  =  —  0  <  4>0  <  n  ,  (5-lc) 

r  r 

where  r  =  y^(I  —  x,<)2  -f-  y2. 

The  field  at  the  receiver  due  to  a  given  ray  path  is  constructed  by  return¬ 
ing  to  the  expression  in  Eq.  (2-16),  rewritten  here  for  convenience: 

,  e,kfll  ik 

p(x,y,z)  =  — —  =  —  /  /  si n0d<j>de  z>0  .  (5-2) 

/ij  Zir  Je= o  J<t>= o 

Equation  (5-2)  expresses  the  field  at  receiver  position  (x.  y,  z)  due  to  a  point  source 
at  the  origin  as  an  integral  over  plane  waves.  Translating  the  source  position 
from  the  origin  to  (x,*,0, z4»),  setting  the  receiver  position  to  (xr,yr,zr),  and 
substituting  for  kx,  kv,  and  k,  from  Eqs.  (2-13),  we  rewrite  Eq.  (5  2)  as 

p(xr,yr,zT)  =  —  /  /  s\n0d<t>d0  ,  (5  3) 

Zn  Je- o  J<t>= o 

where 

$(0,<i>)  =  I-{sin<?((xr  -  x,«)cos<£  +  yr  sin  <£]  -f  (zr  -  zv)cos0)  (5  1) 
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Introducing  <f>0  from  Eq.  (5-lc)  into  the  expression  for  $,  we  obtain 

$($,</>)  =  A:[r  sin0(cos  <^>o  cos  +  sin  <^sin  <p)  +  zcos0]  (5-5) 

=  fc(rsin0cos(</>  —  <f>o)  +  z  cos#]  .  (5-6) 

After  substituting  r  =  R\  sin#o  and  z  =  R\  cos#0,  we  finally  express  4>  as 

^(9,<f>)  =  kRi[s\n0osinOcos((i)  —  <f>0)  +  cos60cos6]  .  (5-7) 

The  next  step  is  to  multiply  the  integrand  of  the  plane  wave  integral 
Eq.  (5  3)  by  the  appropriate  plane  wave  reflection  coefficients.  In  order  to  find 
the  angle  of  incidence  with  which  a  given  plane  wave  strikes  a  given  interface,  we 
use  the  direction  cosine  representation  of  a  vector  v: 

v  :  (cos  [3Z  ,  cos  ,  cos  0ZJ  ,  (5-8) 

where  / 3Z  is  the  angle  between  v  and  the  positive  j-axis,  and  0y  and  0Z  are  defined 
analogously.  The  angle  t/>  between  vectors  vx  and  v2  is  then  given  by 

cos  tp  =  vx  ■  v2  •  (5-9) 

From  simple  geometry,  it  can  be  shown  that  the  direction  cosine  repre¬ 
sentation  of  the  ray  with  angles  6  and  <p ,  as  shown  in  Fig.  5- 1(a),  is 

fii-ay {Q ,  <t>) '■  (sin  0  cos  <f> ,  sin  0  sin  <f> ,  cos  O^j  .  (5-10) 

This  ray  is  perpendicular  to  the  wave  fronts  of  the  associated  plane  wave.  We  also 
require  the  representation  of  the  normals  to  the  interfaces.  For  the  ;th  bottom 
reflection  (as  counted  from  the  receiver),  the  normal  to  the  imaged  ocean  bottom 
has  direction  cosines 

Vj  :  (sin [(2.7  -  l)a]  ,  0  ,  ±cos[(2y  -  1  )o]j  (5-11) 

The  positive  sign  is  taken  for  imaged  bottoms  that  lie  below  the  z  =  0  plane- 
frays  traveling  upward  at  the  receiver),  and  the  negative  sign  is  taken  for  imaged 
bottoms  that  lie  above  the  c  =  0  plane  (rays  traveling  downward  at  the  receiver). 
The  signs  are  taken  so  that  the  normal  to  each  interface  points  in  the  same  direc¬ 
tion  as  the  passage  of  the  ray  from  source  to  receiver  (~jr/2  <  t/’j  <  7r/2).  From 
Eq.  (5  9),  the  incident  angle  i/'y  of  a  plane  wave,  with  angles  0  and  <p,  at  the  yth 
bottom  reflection  is  then  given  by 

cos  —  sin  0  cos  0  sm((2y  —  1  )o]  i  cos  0  cos[(2/  —  1  )<  y  ]  (5  12) 
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Note  that  tpj  =  0  ±  (2 j  —  1)q  for  the  two-dimensional  wedge  problem,  where 
cos  <t>  —  ±1. 

We  now  construct  the  integral  representation  of  the  ray  field  in  the  three- 
dimensional  wedge  for  the  ray  path  with  Nt  surface  reflections  and  Nf,  bottom 
reflections.  To  do  so,  we  multiply  the  integrand  of  Eq.  (5-3)  by  the  reflection 
coefficients  V  evaluated  at  the  angles  given  in  Eq.  (5-12): 

p(«o,  *>,*»)  =  (-1)"*^  t_ 0  Z*™  (n  sin  0 d<f> d0  , 

(5-13) 

where  $  is  given  in  Eq.  (5-7).  Unlike  the  two-dimensional  case,  we  cannot  per¬ 
form  the  integration  over  <j>  analytically  because  the  incident  angles  ipj  are  now 
dependent  on  <t>  as  well  as  0.  This  dependence  is  an  indication  that  the  spreading 
factors  of  the  ray  tube  in  the  9-  and  ^-directions  are  not  independent. 

As  we  have  done  with  the  previous  plane  wave  integrals,  we  move  the 
reflection  coefficient  terms  into  the  exponential  factor  and  write 

p(0o,4o,Ri)  =  (-l)N‘ £  /~,0°  P  e'HW  sine  dtdO  ,  (5-14) 

Zn  j$= o  J<t>=o 

where 

p  =  kRx  (5-15) 

/ ( 6 ,  <i>)  =  i[sin  0O  sin  6  cos (<f>  —  <j> o)  +  cos  60  cos  6} 

+  !'[««.«]  •  (5-16) 

Note  that  the  phase  function  /  is  now  a  function  of  two  angular  variables  9  and  <f>. 
Application  of  the  method  of  steepest  descent  in  this  case  is  not  straightforward 
because  (1)  there  are  two  variables  of  integration  and  (2)  all  of  the  limits  of 
integration  are  not  infinite. 


5-3  AN  APPROXIMATE  METHOD  OF  SOLUTION 


We  now  propose  an  approximate  method  for  evaluating  the  ray  field  given 
in  Eq.  (5-14).  The  first  step  is  to  find  all  the  saddle  points  of  the  function  /  in 
Eq.  (5-16),  where  a  saddle  point  now  consists  of  a  pair  of  possibly  complex  angles 
defined  by 


gWf)  _  df[0,<t>) 

89  d<f> 


at  9  =  7  ,  <{>  =  x- 


(5-17) 
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In  order  to  be  able  to  perform  the  integration  over  <p  for  each  saddle  point,  we 
evaluate  the  reflection  coefficient  terms  at  <f>  =  \  and  take  them  outside  the 
integral  over  <p.  This  approximation  yields  the  integral 


p{0o,<f>o,  Ri)  = 


y*  ,0°  jy  *  eikR> *»»«o«n0co«(*-$o) 

exp^/?!  ^icos^ocos^  -f  —■  £  In  V\tpj(0, \)]^ 


si  n0d9  , 

(5-18) 


and  allows  the  integral  over  <f>  to  be  performed  analytically  as  in  Eq.  (2  22).  At 
this  point,  the  remainder  of  the  derivation  given  in  Sec.  2  2  follows.  The  result 
is  an  expression  very  similar  to  that  constructed  for  the  two-dimensional  wedge 
problem  in  Eq  (3-9),  except  that  the  reflection  coefficients  are  evaluated  at  angles 
ipj,  which  depend  on  the  saddle  point  angle 


p{0o,  Ri) 


This  integral  may  be  solved  using  the  same  saddle  point  methods  used  in  Sec.  3-  5 
for  the  two-dimensional  wedge. 

A  problem  arises,  however,  when  applying  the  saddle  point  method  to 
Eq.  (5-19),  where  the  function  /  has  the  form 


i  Nb 

f{O)  =  icos(0-Oo)  +  —  X>V(<M0iX)]  -  (5-20) 

i  ]=i 


Note  that  the  saddle  point  7  found  from  the  function  /  in  Eq.  (5-16)  is  a  saddle 
point  for  /  in  Eq.  (5-20)  only  if  \  =  <po ■  This  condition  may  not  be  satisfied 
exactly,  especially  if  the  eigenray  undergoes  displacements  at  the  bottom.  We 
could  find  a  new  saddle  point  7  for  /,  or  we  could  keep  the  old  saddle  point  from 
/  by  replacing  /  in  Eq.  (5-20)  by 


f{0)  =  f{0,  <p  -  x)  =  i[sin0osin0cos(\  -  d>u)  +  cos0ocos0] 
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The  latter  choice  appears  to  make  more  sense  since  the  saddle  points  7  and  \  that 
satisfy  the  original  saddle  point  criterion  are  preserved. 

Some  of  the  details  of  the  calculations  required  for  this  problem  will  now 
be  given.  To  find  the  saddle  points  defined  in  Eq.  (5-17),  the  partial  derivatives 
of  /  in  Eq.  (5-16)  must  be  calculated.  The  required  expression  are 


df 

30 


dj 

3<f> 


t[sin  0q  cos  0  cos(<j)  —  <j> 0)  —  cos  0O  sin  6) 

.......  ,  .  ,  1  £  v\u»A)\Wi 


(5-22) 

(5-23) 


where  the  partial  derivatives  of  ipj ,  obtained  from  Eq.  (5-12),  are: 

=  —  (cos#cos<^sin[(2_;  —  1  )a]  T  sin  0  cos[(2j  —  1  )a])  /  sini/>;  (5-24) 

30  '  ' 

=  (sin^sin^sin[(2j  -  ljaj^/sint/’j  .  (5-25) 

The  algorithm  for  finding  saddle  points  would  use  Taylor  series  expansions  of  a 
function  of  two  variables  rather  than  one. 


5-3.1  An  analysis  of  the  approximation 

By  removing  V[rpj(0,<f>)\  from  the  integral  over  <f>  in  Eq.  (5-18),  we  are 
effectively  ignoring  any  coupling  between  ray-tube  spreading  factors  in  the  0-  and 
^-directions.  In  reality,  this  coupling  does  occur  for  rays  that  do  not  travel  straight 
upslope  or  downslope  (those  not  confined  to  the  y  —  0  plane  in  Fig.  5-1 ).  Coupling 
between  ray-tube  spreading  factors  corresponds  to  the  phenomenon  known  as 
horizontal  refraction,  where  the  path  of  a  ray  as  viewed  from  above  the  ocean 
surface  appears  to  bend  as  the  ray  undergoes  successive  reflections.  Horizontal 
refraction  occurs  when  a  ray’s  plane  of  incidence  (the  plane  containing  the  ray 
and  the  normal  to  the  interface)  changes  at  successive  bottom  interactions.  The 
approximation  we  have  introduced  may  result  in  errors  in  the  calculation  of  the 
two  eigenray  spreading  factors  G$  and  G $  defined  in  Sec.  2-7.1,  but  the  eigenray 
paths  themselves  are  not  affected  since  the  full,  two-dimensional  saddle  point 
criterion  in  Eq.  (5-17)  is  used. 

The  validity  of  the  approximation  made  to  arrive  at  Eq.  (5-18)  depends 
on  how  strongly  V  varies  with  <f>  in  the  neighborhood  of  the  saddle  points  7  and 
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Since  V'  is  large  in  the  neighborhood  of  the  critical  angle,  the  approximation 
is  potentially  poorest  when  an  eigenray  undergoes  a  lateral  wave  reflection  on  one 
of  its  bottom  interactions  (where  ^  as  6a).  But  since  the  lateral  wave  rays  make 
up  a  small  part  of  the  total  field,  slight  errors  in  their  evaluation  should  not  harm 
the  total  calculation  significantly.  It  should  also  be  noted  that  when  propagation 
is  mainly  upslope  or  downslope,  V’j  does  not  vary  significantly  with  <j>  (as  will 
be  shown  in  the  next  paragraph),  so  the  lateral  wave  calculations  should  exhibit 
smaii  errors  for  these  cases. 

For  the  more  important  rays  that  have  saddle  points  somewhat  removed 
from  0CT,  the  validity  of  the  approximation  introduced  in  Eq.  (5-18)  depends  on 
how  strongly  the  incident  angles  ipj  vary  with  <j>.  Thus,  the  magnitude  of  dxl>} / ()<t> 
from  Eq.  (5-25)  should  be  small.  We  see  that  this  is  true  when’  (1)  7  is  small, 
(2)  X  is  small,  or  (3)  (2 j  —  1)q  is  small.  Again,  the  second  condition  is  satisfied 
for  mainly  upslope  or  downslope  propagation.  The  third  condition,  which  holds 
for  small  wedge  angles  and  few  bottom  reflections,  can  be  satisfied  even  in  the 
cross-slope  case;  for  example,  if  $a  is  not  small  and  the  range  between  source  and 
receiver  is  not  too  large,  then  eigenrays  with  large  Nb  are  not  important  because 
they  suffer  multiple  partial  reflections. 

In  cases  where  eigenrays  undergo  significant  horizontal  refraction,  the 
approximation  we  have  outlined  here  may  not  be  acceptable.  In  such  cases,  it 
appears  that  a  two-dimensional  integration  must  be  performed.  One  problem 
that  arises  when  attempting  to  formulate  the  integrals  in  terms  of  the  saddle  point 
method  is  that  the  limits  of  integration  in  Eq.  (5-14)  are  not  infinite.  When  the 
integral  over  <j>  cannot  be  performed  analytically,  it  appears  difficult  to  manipulate 
the  integral  over  6  such  that  the  contour  starts  at  —  tt/2  +  too  instead  of  0.  It 
remains  to  be  investigated  if  the  method  of  steepest  descent  can  be  applied  to  the 
integral  whose  contour  starts  at  0. 


'The  denominator  of  Eq.  (5-25)  can  make  large  only  if  is  dose  to  zero.  Hut  this 

condition  implies  that  the  eigenray  must  be  traveling  straight  upslope  or  downslope,  in  which 
case  we  have  shown  that  the  approximation  is  valid. 


CHAPTER  6 


THE  MODELING  OF  A  DIRECTIONAL  SOURCE 


6-1  INTRODUCTION 

The  methods  described  herein  for  modeling  a  uniform  acoustic  point 
source  may  be  extended  to  a  directional  source  by  allowing  the  spatial  coordi¬ 
nates  of  the  point  source  to  take  on  complex  values.  It  was  first  observed  by 
Deschamps45  that  the  solution  exp (ikR\)/ R\  to  the  Helmholtz  equation  represents 
the  field  of  a  Gaussian  beam  when  the  distance  Rt  is  calculated  from  a  complex 
source  position  to  a  real  receiver  position.  By  analytic  continuation,  known  exact 
or  asymptotic  solutions  to  point  source  problems  can  then  be  extended  to  direc¬ 
tional  source  problems  by  simply  introducing  the  appropriate  complex  quantities. 
This  method  is  used  in  Refs.  46  and  47,  where  asymptotic  forms  for  the  reflection 
and  transmission  of  electromagnetic  beams  at  a  dielectric  interface  are  given  and 
discussed. 

6-2  THE  FIELD  DUE  TO  A  COMPLEX  POINT  SOURCE  IN  FREE 
SPACE 

We  begin  by  describing  the  nature  of  the  field  produced  by  a  point  source 
in  complex  space.  Referring  to  Fig.  6-1,  we  initially  assume  the  point  source  to 
be  located  at  Cartesian  coordinates  (0, 0,  *6)  and  consider  the  field  at  coordinates 
(x,  0,  z).  (We  have  restricted  our  attention  to  the  x-z  plane  with  the  understanding 
that  the  field  is  cylindrically  symmetric  about  the  2-axis.)  The  non-zero  imaginary 
part  6  of  the  source  coordinate  results  in  a  beam  whose  axis  is  in  the  2-direction. 
Denoting  the  source-receiver  distance  /?],  we  find 

Hi  -  \J{z  ~  ib)7  +  J2  ,  (6-1) 

where  R\  is  now  complex.  Since  we  are  most  often  interested  in  the  field  near  the 
axis  of  a  directional  source,  we  concentrate  on  the  paraxial  region  «  ;2  +  b2. 
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Figure  6-1  Geometry  for  the  directional  source.  A  point  source  is  placed  at 
(*,y,z)  COORDINATES  (0, 0,  *6) ,  WHICH  PRODUCES  A  BEAM  IN  THE  POSITIVE  2-DIRECTION. 


where  we  can  make  the  approximation 


Ri 


z  -  ib  + 
z  —  ib  + 


2 (z  —  ib) 
x2(z  +  ib) 
2{z2  +  b>) 


(6-2) 


Using  this  expression  for  ft,  in  the  exponential  of  the  field  due  to  a  point  source, 
we  find 


e'kR'  ss  exp  (kb)  exp  ( ikz  1  + 


jtfcz  1 


2(z2  +  b2) 


]} 


exp 


kb 


2  (z2  +  62). 


(6-3) 


The  first  exponential  factor  on  the  right  hand  side  of  Eq.  (6-3)  is  a  constant  that 
can  be  very  large  or  small,  so  we  normalize  the  field  by  the  factor  exp(-kb).  The 
second  exponential  factor  gives  the  mathematical  form  for  the  phase  fronts  in  the 
paraxial  region.  The  final  exponential  factor  represents  a  Gaussian  amplitude 
distribution  of  the  beam  about  the  z-axis. 

The  square  root  operation  required  to  calculate  the  complex  distance  ft, 
in  Eq.  (6-1)  introduces  a  branch  cut  in  real  space.  If  we  do  not  limit  ourselves  to 
the  x  —  z  plane  and  define  the  radicand  as 


s(x,y,z)  =  (z  -  ib)2  +  x2  +  y2  =  (z2  -  62  +  x2  +  y2)  -  i(26z)  ,  (6  1) 

then  the  branch  points  are  located  at  (x,y,z)  values  where  both  the  real  and 
imaginary  parts  of  s  vanish.  It  is  easily  seen  that  the  branch  points  form  a  circle 
of  radius  6  in  the  z  =  0  plane  defined  by  x2  +  y2  =  62.  We  define  the  branch 
cuts  for  the  square  root  operation  to  be  at  points  where  Im[s]  =  0  and  Refs]  <  0. 
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The  first  condition  implies  that  z  —  0,  and  the  second  condition  implies  that 
x2  +  y2  <  b7.  The  branch  cut  surface  is,  therefore,  a  disk  of  radius  b  in  the  x  —  y 
plane,  lying  in  real  space  perpendicular  to  the  beam  axis  (see  Fig.  6-1). 

With  the  branch  cut  as  defined  above,  the  evaluation  of  /?,  can  be  made 
either  on  the  upper  Riemann  sheet,  for  which  Re[/2j ]  >  0,  or  on  the  lower  Riemann 
sheet,  for  which  Re[R,]  <  0.  When  we  include  the  suppressed  time  dependence 
exp (—iu>t),  we  see  that 

— e'(kWi_w<)  =  -  c-»m[flilc«(fcRe[R,l-^)  (6-5) 

R\  R\ 

From  the  last  exponential  factor  in  Eq.  (6-5),  we  see  that  on  the  upper  Riemann 
sheet  we  have  an  outgoing  wave,  while  on  the  lower  Riemann  sheet  we  have  an 
incoming  wave.  The  significance  of  the  branch  cut  is  that  the  field  is  discontinuous 
across  the  cut  if  all  calculations  of  /?i  are  made  on  the  same  Riemann  sheet. 
Therefore,  points  on  the  branch  cut  disk  in  real  space  may  be  viewed  as  equivalent 
sources  that  together  represent  the  single  point  source  in  imaginary  space. 


6-3  DETERMINATION  OF  COMPLEX  SOURCE  PARAMETERS 

The  method  described  here  for  modeling  a  directional  source  may  be  used 
to  investigate  the  behavior  of  beams  from  a  theoretical  standpoint  or  to  model 
an  actual  acoustic  source.  In  determining  the  parameters  of  the  complex  source 
to  be  used,  it  is  appropriate  to  consider  two  basic  characteristics  of  the  beam  to 
be  modeled:  (1)  a  measure  of  the  beamwidth  and  (2)  a  measure  of  the  phase 
front  curvature  of  the  beam  along  its  axis.  From  these  two  quantities  we  derive 
expressions  for  z  and  b.  The  latter  can  be  thought  of  as  the  distance  into  imaginary 
space  the  source  is  placed,  while  the  former  is  the  axial  distance  in  real  space  from 
the  source  to  the  region  where  the  beam  has  the  desired  characteristics. 

We  define  the  spatial  “variance”  or  width  of  the  Gaussian  beam  described 
in  the  last  exponential  factor  of  Eq.  (6-3)  by 


For  x  =  ±cr  the  field  is  down  by  e  ‘/2  =  0.6065,  or  —4.34  dB,  from  its  peak  at 
i  =  0;  for  x  =  ±2<t  the  field  is  down  by  —17.32  dB. 

The  radius  of  curvature  of  the  phase  fronts  in  the  paraxial  region  can 
be  found  from  the  second  exponential  factor  in  Eq.  (6  3).  The  equation  for  the 


curves  of  constant  phase  in  the  x-z  plane  is 


K1  +  2F;i,Tj)=const  •  (6“7) 

The  radius  of  curvature  of  these  phase  fronts  on  the  z-axis  may  be  found  by 
considering  a  point  on  the  front  a  small  distance  Ax  from  the  z-axis,  drawing 
the  tangent  to  the  wavefront  at  that  point,  and  finding  the  intersection  of  the 
perpendicular  to  the  tangent  with  the  z-axis  (see  Fig.  6-2).  The  distance  from 


Figure  6-2  Calculation  of  the  radius  of  phase  front  curvature  Rc 


the  intersection  point  to  the  phase  front  in  the  limit  of  small  Ax  is  then  the  radius 
of  curvature  Rc.  Some  simple  calculations  based  on  Fig.  6-2  reveal  that 


Rc  z"(x  =  0) 


(6  8) 


where  z(x)  is  the  functional  form  of  the  phase  front  and  z"(x  =  0)  is  its  second 
derivative  with  respect  to  x  evaluated  at  x  =  0.  After  taking  (implicitly)  two 
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derivatives  of  Eq.  (6-7)  with  respect  to  x  and  evaluating  the  resulting  function  at 
x  =  0,  we  find  that  the  radius  of  curvature  of  the  phase  fronts  on  the  z-axis  is 


Rc  = 


z3  +  b3 


(6-9) 


Equations  6-6  and  6-9  can  be  solved  simultaneously  for  the  two  variables 

z  and  b: 


(ko2)2  Rc 
R2  +  (ka*)* 

{ka2)R2 
R2  +  (ka2y 


(6-10) 

(6-11) 


Thus,  given  the  quantities  a  and  Rc  that  characterize  a  beam  in  a  certain  region 
of  interest,  one  can  solve  for  z  and  b.  To  model  the  beam  by  a  complex-valued 
point  source,  one  would  place  the  point  source  a  distance  z  in  real  space  and  a 
distance  b  in  imaginary  space  away  from  the  region  of  interest. 


6-4  REFLECTED  AND  TRANSMITTED  FIELDS  DUE  TO  AN  IN¬ 
CIDENT  BEAM 

The  problem  in  which  we  are  interested  is  that  of  modeling  the  reflection 
and  transmission  of  an  acoustic  beam  at  a  fluid-fluid  interface  lying  in  the  z  =  0 
plane.  We  consider  a  beam  that  strikes  the  interface  at  the  origin  of  our  coordinate 
system  with  angle  of  incidence  9g  (see  Fig  6-3)  Values  of  the  h*»»mwidt,h  a  and 
the  radius  of  phase  front  curvature  Rc  in  the  region  of  the  origin  are  to  be  specified 
or  estimated  from  experimental  measurements.  Equations  (6-10)  and  (6-11)  are 
then  used  to  calculate  z  (which,  to  avoid  confusion  with  the  z-axis  in  the  present 
case,  we  will  refer  to  as  p0 )  and  6.  After  the  appropriate  rotation  and  translation 
of  the  coordinate  system,  the  point  source  is  placed  at  [x,z]  coordinates 

(-po  +  i&)sin0B,(-po  +  i6)cos0B]  •  (6-12) 

Note  that  we  cannot  allow  the  disk  of  radius  b  perpendicular  to  the  beam  axis 
to  intersect  the  interface  or  else  the  equivalent  sources  would  extend  into  the 
bottom.  This  restriction  requires  bsinOg  <  p0cos9g.  The  result  is  that  it  may 
not  be  possible  to  model  a  narrow  beam  (6  large)  incident  at  small  grazing  angles 
( 0g  near  90°). 
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Figure  6-3  Geometry  for  the  reflection  and  transmission  of  a  beam  produced 

BY  A  POINT  SOURCE  IN  COMPLEX  SPACE. 
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The  method  for  finding  the  reflected  and  transmitted  fields  due  to  the 
directional  source  is  the  same  as  for  the  point  source,  except  that  certain  geometri 
cal  quantities  are  complex.  Since  the  source  coordinates  have  been  made  complex, 
the  distances  r  and  Rx  and  the  incident  angle  0O  in  the  reflected  field  expression 
Eq.  (2-46)  and  in  the  transmitted  field  expression  Eq.  (2  75)  are  also  complex. 
The  fact  that  Rx  is  now  complex  requires  a  slight  redefinition  of  the  quantities  p 
and  f(9)  associated  with  the  method  of  steepest  descent  in  Kq.  (2-29).  For  the 
directional  source  case  we  call  these  quantities  p  and  f(9)  .  In  order  to  keep  /> 
real,  we  set 

p=\p\  =  \kRi\  (6  id) 

and 

m  =  -f(0)  =  ^no)  .  (fin) 

p  i*«ii 

where  f(9)  is  defined  in  Eq.  (2  47b)  for  the  reflection  case  and  in  Kq.  (2  76)  for 
the  transmission  case.  Note  that  Eqs.  (613)  and  (6  14)  do  not  constitute  an 
approximation  because  pf  —  pf.  As  before,  we  find  the  saddle  point(s)  y  as  the 
valuc(s)  of  9  for  which  f'{9)  vanishes.  The  field  is  then  evaluated  by  either  the 
saddle  point  approximation  or  by  numerical  integration  along  the  path  of  steepest 
descent. 

6-5  EXAMPLES  OF  COMPUTED  FIELDS 

Examples  of  the  acoustic  field  due  to  the  reflection  and  transmission  of 
a  beam  incident  at  a  fluid-fluid  interface  are  shown  in  Figs.  6  4  to  6  6.  The 
figures  show  contour  lines  of  constant  pressure  amplitude  of  the  reflected  field 
(2  <  0)  and  the  transmitted  field  (:  >  0)  for  beams  incident  on  the  fluid-fluid 
interface  at  angles  9 «  =  55°,  62°.  and  68°.  The  beam  axis  intersects  the  interface 
[characterized  by  n  =  0.866  ( 9ct  =  60°)  and  m  =  2.0]  at  the  origin.  The  beam 
is  specified  to  have  a  variance  of  0  =  5A  and  a  radius  of  wavefront  curvature  of 
Rc  =  50A.  These  characteristics  result  an  imaginary  source  coordinate  of  b  = 
14. 5A  and  p0  =  45.4 A.  Note  that  the  fields  across  the  interface  in  Figs.  6  4  to  6  6 
are  discontinuous.  This  is  because  the  incident  field  due  to  the  beam  source  has 
not  been  added  to  the  reflected  field  above  the'  interface.  To  'est  the  computer 
algorithm  that  produced  the  contour  plots,  the  total  field  above  the  interface  was 
calculated,  and  there  was  less  than  a  0.01  dB  discontinuity  across  the  interface. 
The  condit  ion  of  pressure  cont  inuity  is  one  of  the  two  boundary  conditions  imposed 
at  a  fluid  fluid  interface. 
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The  fields  shown  in  Figs.  6-4  to  6-6  are  intended  to  illustrate  some  ba¬ 
sic  characteristics  of  beam  interaction  at  a  fluid-fluid  interface  when  the  angle  of 
incidence  is  near  the  critical  angle.  Rather  than  attempt  to  undertake  a  compre¬ 
hensive  study  of  this  phenomenon,  we  mean  to  show  that  (1)  the  complex  source 
approach  can  be  used  with  few  restrictions  to  find  the  field  due  to  a  directional 
source,  (2)  basic  characteristics  of  the  reflected  and  transmitted  fields,  such  as 
the  beam  angle  and  beam  displacement  at  the  interface,  can  be  analyzed,  and 
(3)  the  relative  importance  of  different  contributions  to  the  total  field  at  individ¬ 
ual  points  can  be  evaluated.  With  regard  to  the  second  item,  the  axes  of  the 
reflected  and  transmitted  beams  have  been  estimated  by  roughly  connecting  the 
extremal  points  on  successive  contour  lines.  The  beam  displacements,  which  we 
call  A  for  reflection  and  A(  for  transmission,  are  found  by  extending  the  beam 
axes  back  to  the  interface.  Although  this  procedure  is  not  rigorous,  it  does  give  a 
good  estimate  of  the  beam  characteristics. 

In  Fig.  6-4  the  beam  is  incident  on  the  interface  at  an  angle  of  55°,  which 
is  5°  steeper  (smaller)  than  critical  incidence.  In  this  case  both  the  reflected  and 
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Figure  6  4  Reflected  and  transmitted  field  for  an  incident  beam  with  ancle 
9b  =  55° 

transmitted  beams  are  well  defined,  even  near  the  interface.  The  reflected  beam 
has  an  apparent  displacement  of  A  =  2.6A,  while  the  transmitted  beam  has  an 
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apparent  displacement  of  At  =  2.0A.  These  results  may  seem  surprising  at  first 
because  the  formulas  in  Eqs.  (2-60)  and  (2-81)  indicate  that  the  ray  displacements 
at  these  angles  should  be  zero.  However,  it  appears  that  these  phenomena  may  be 
explained  consistently  by  considering  the  facts  that  the  incident  beam  is  actually 
composed  of  a  spectrum  of  plane  waves  about  the  beam  axis  and  that  rays  drawn 
perpendicular  to  these  plane  waves  strike  the  interface  at  different  points.  For 
example,  the  rays  incident  at  0  =  55°  strike  the  interface  at  the  origin  and  are 
reflected  with  a  plane  wave  reflection  coefficient  V  of  0.61,  while  those  at  the 
critical  angle  9  —  60°  strike  the  interface  to  the  right  of  the  origin  and  are  reflected 
with  a  coefficient  of  1.0.  Although  the  latter  rays  are  initially  weaker  than  the 
former  due  to  the  fact  that  they  are  5°  off  the  beam  axis,  their  larger  reflection 
coefficient  could  make  them  stronger  after  reflection.  This  may  explain  not  only 
why  the  reflected  beam  axis  has  been  shifted  to  the  right,  but  also  why  its  angle 
is  slightly  larger  (57°)  than  that  of  the  incident  beam. 

In  the  case  of  the  transmitted  beam  in  Fig.  6-4,  one  must  also  consider  the 
spreading  introduced  by  refraction  according  to  Snell’s  law.  For  incident  angles 
slightly  steeper  than  critical,  a  small  change  in  incident  angle  results  in  a  large 
change  in  transmitted  angle.  The  resulting  spreading  is  largest  near  the  critical 
angle  and  diminishes  toward  steeper  angles.  This  factor  may  explain  why  the 
transmitted  beam  in  Fig.  6-4  has  a  measured  angle  of  67.4°,  while  Snell’s  law 
predicts  an  angle  of  71.1°  for  an  incident  beam  at  55°.  The  shift  to  a  steeper 
angle  may  be  due  to  the  decreased  spreading  upon  refraction  at  steeper  angles. 
However,  since  the  steeper  rays  strike  the  interface  to  the  left  of  the  beam  axis, 
one  might  expect  the  displacement  upon  transmission  to  be  negative  rather  than 
the  indicated  positive  value.  Clearly,  the  present  arguments  are  inadequate  to 
account  for  both  the  displacement  and  angle  of  the  transmitted  beam. 

Figure  6-5  shows  the  reflected  and  transmitted  fields  for  a  beam  incident 
at  62°,  which  is  2°  larger  than  the  critical  angle.  The  reflected  beam  is  displaced 
by  4.9A  and  has  an  angle  slightly  larger  than  the  incident  beam.  These  phenomena 
can  be  explained  using  the  same  argument  as  for  the  55°  beam:  the  effective  beam 
center  is  shifted  to  shallower  angles  and  to  the  right  because  the  steeper  angles 
experience  more  loss  upon  reflection.  The  ray  displacement  at  62°,  according  to 
Eq.  (2-60),  should  be  A  =  1.9A.  The  measured  value  of  A  =  4.9A  is  larger  due 
to  the  additional  effect  of  the  beam  center  being  shifted  to  the  right. 

The  transmitted  beam  in  Fig.  6-5  is  displaced  1 .8A  to  the  left  and  has 
an  r  >gle  of  75.4°.  By  checking  the  eigenray  list  for  receiver  positions  along  the 
beam  axis,  it  is  apparent  that  the  main  contribution  to  the  field  is  from  eigenrays 
that  are  incident  on  the  interface  at  angles  slightly  steeper  than  the  critical  angle. 
These  eigenrays  represent  a  part  of  the  beam  that  strikes  the  interface  to  the  left 
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Figure  6-5  Reflected  AND  transmitted  field  for  an  incident  beam  with  ancle 
9  b  =62°. 

of  the  main  beam  axis.  Therefore,  the  apparent  displacement  of  the  transmitted 
beam  is  negative.  The  beam  angle  in  the  bottom  is  due  to  a  combination  of  the 
following  factors:  (1)  the  decreasing  spreading  loss  upon  refraction  as  the  incident 
angle  becomes  steeper  than  critical,  (2)  the  decreasing  strength  of  the  incident 
beam  for  angles  steeper  than  the  beam  axis,  and  (3)  the  decreasing  transmission 
coefficient  [see  Eq.  (2-73)]  as  the  incident  angle  becomes  steeper  than  critical.  The 
transmitted  beam  at  62°  incidence  is  shallower  than  at  55°  incidence,  as  expected. 

The  final  example  appears  in  Fig.  6-6,  where  the  incident  beam  angle 
is  68°.  The  reflected  beam  has  an  angle  of  65.6°,  compared  to  the  value  of  68° 
predicted  from  Snell’s  law,  and  a  measured  displacement  of  1.1  A,  compared  to  a 
predicted  ray  displacement  of  1.8A.  By  virtue  of  the  fact  that  the  incident  beam  is 
8°  shallower  than  critical,  one  would  expect  there  to  be  little  shift  of  the  effective 
beam  center  due  to  unequal  reflection  on  opposite  sides  of  the  incident  beam  axis. 
One  possible  explanation  is  that  the  deviation  is  due  to  interference  from  the 
lateral  wave  field.  A  check  of  the  eigenray  list  shows  that  the  lateral  wave  ray  is 
typically  30  dB  down  from  the  totally  reflected  ray  along  the  beam  axis.  On  the 
interface  side  of  the  beam  axis  the  lateral  wave  ray  becomes  relatively  stronger. 
This  may  account  for  the  apparent  asymmetry  of  the  reflected  beam. 

The  transmitted  beam  in  Fig.  6  6  has  an  angle  of  80.5°  and  a  negative 
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Figure  6  6  Reflected  and  transmitted  field  for  an  incident  ream  with  angle 
9b  =68°. 


displacement  of  14. 2A.  The  displacement  to  the  left  is  larger  in  this  case  than 
the  previous  because  energy  incident  at  angles  steeper  than  critical  strikes  the 
interface  farther  to  the  left.  The  transmitted  beam  is  significantly  weaker  in  this 
case  because  the  penetrating  energy  is  farther  off  the  beam  axis.  As  expected, 
incident  beams  at  shallower  angles  excite  transmitted  beams  at  shallower  angles. 

In  both  Figs.  6-5  and  6-6,  the  evanescent  field  in  the  upper  wavelength 
or  so  of  the  bottom  is  evident.  In  the  latter  case  an  interference  pattern  between 
this  field  and  the  transmitted  field  can  be  seen,  particularly  at  a  range  of  4A  and  a 
depth  of  1A.  These  effects  can  make  it  impossible  to  trace  a  beam  axis  all  the  way 
back  to  the  interface.  The  transmitted  beam  in  such  cases  is  only  distinguishable 
a  sufficient  distance  from  the  interface. 

6-6  SADDLE  POINT  STRUCTURE 

The  saddle  points  and  their  associated  steepest  descent  paths  for  receivers 
along  several  horizontal  lines  are  shown  in  Figs.  6-7  and  6-8.  In  each  figure  both 
the  reflected  and  transmitted  fields  two  wavelengths  from  the  interface  (rr  =  ±2A) 
are  sampled  ten  times  from  horizontal  ranges  of  — 5A  to  25A.  Figure  6  7  is  for  the 
incident  beam  at  55°  shown  in  Fig.  6  4,  and  Fig.  6  8  is  for  the  incident  beam  at 
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68°  shown  in  Fig.  6-6. 

The  general  effect  of  the  source  directionality  on  the  saddle  point  loca¬ 
tions  is  to  shift  the  saddle  points  farther  off  the  real  0-axis  as  the  eigenray  angle 
moves  farther  off  the  beam  axis.  The  computer  algorithm  for  finding  the  cor¬ 
rect  saddle  points  must  be  more  sophisticated  in  the  present  case  because  there 
are  more  possibilities  for  how  the  saddle  points  migrate  and  evolve  and  how  the 
correct  contour  of  integration  is  completed.  For  example,  in  Fig.  6-7(a)  the  par¬ 
tially  reflected  eigenrays  evolve  smoothly  into  the  lateral  wave  eigenrays,  whereas 
they  usually  evolve  into  the  totally  reflected  eigenrays,  as  they  do  in  Fig.  6-8(a). 
Another  possible  scenario  is  that  for  the  transmitted  field,  the  steepest  descent 
path  of  the  transmitted  ray  can  loop  around  the  critical  angle,  but  both  endpoints 
terminate  towards  zr/2  —  zoo  instead  of  —  tt/2  +  zoo.  In  this  case,  an  evanescent 
ray  saddle  point  must  be  found  on  the  lower  Riemann  sheet.  Its  steepest  descent 
path  completes  the  contour  of  integration. 

In  Fig.  6-7(a)  we  see  that  the  reflected  field  at  the  largest  three  ranges 
requires  two  saddle  points,  while  the  other  ranges  require  only  one.  The  lateral 
wave  rays  are  about  12  dB  weaker  than  the  totally  reflected  rays.  In  part  (b) 
of  the  figure  we  see  that  the  transmitted  field  at  the  largest  two  ranges  requires 
two  saddle  points:  one  transmitted  and  one  evanescent.  The  latter  rays  are  about 
40  dB  weaker  than  the  former. 

In  Fig.  6-8(a)  for  the  reflected  field,  we  see  that  all  but  the  first  receiver 
have  a  lateral  wave  saddle  point.  Note  how  the  saddle  points  cross  the  real  0-axis 
around  68°,  which  is  the  angle  of  the  incident  beam.  For  the  transmitted  field  in 
part  (b)  we  see  that  all  but  the  first  receiver  have  the  evanescent  ray  as  well  as 
the  transmitted  ray.  The  evanescent  field,  however,  is  20-29  dB  weaker  than  the 
transmitted  field,  despite  the  fact  that  the  transmitted  rays  near  60°  are  quite  far 
off  the  incident  beam  axis. 

6-7  CONCLUSIONS 

In  this  chapter  we  have  demonstrated  how  our  method  for  modeling 
acoustic  propagation  in  terms  of  eigenrays  can  be  modified  to  accommodate  direc¬ 
tional  sources.  We  have  described  the  nature  of  the  directional  field  produced  by 
placing  a  point  source  in  complex  space,  and  we  interpret  the  associated  branch 
cut  surface  as  a  disk  of  sources  in  real  space  equivalent  to  the  point  source  in 
complex  space.  One  method  for  modeling  a  beam  with  a  given  radius  of  phase 
.ront  curvature  and  beamwidth  has  been  proposed.  Finally,  we  have  given  sev¬ 
eral  examples  of  the  reflected  and  transmitted  fields  due  to  a  beam  incident  on  a 
fluid-fluid  interface  near  the  critical  angle.  The  same  techniques  described  in  this 


104 


Figure  6-7  Saddle  points  and  steepest  descent  paths  for  (a)  the  reflected  field 
2A  ABOVE  THE  INTERFACE  AND  (b)  THE  TRANSMITTED  FIELD  2A  BELOW  THE  INTERFACE  FOR 
THE  55°  INCIDENT  BEAM  OF  FlC.  6  4. 
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Figure  6-8  Saddle  points  and  steepest  descent  paths  for  (a)  the  reflected  field 
2A  ABOVE  THE  INTERFACE  AND  (b)  THE  TRANSMITTED  FIELD  2A  RELOW  THE  INTERFACE  FOR 
THE  68°  INCIDENT  BEAM  OF  FlC.  6  6. 
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chapter  for  modeling  a  beam  incident  on  a  single  interface  could  be  used  to  model 
a  directional  source  in  a  flat  or  sloping  waveguide  by  simply  using  the  appropriate 
plane  wave  integrals  given  in  Chapter  3. 


CHAPTER  7 


CONCLUSIONS 


The  contribution  of  the  research  presented  in  this  dissertation  is  the  im¬ 
provement  of  ray  techniques  to  model  acoustic  propagation  in  shallow  water.  In 
Chapter  2  we  developed  a  method  for  computing  the  Held  due  to  a  point  source 
in  the  presence  of  a  fluid-fluid  interface.  After  formulating  plane  wave  integral 
representations  of  the  reflected  and  transmitted  fields,  we  use  a  variation  on  the 
method  of  steepest  descent  to  solve  the  integrals.  By  allowing  the  reflection  and 
transmission  coefficients  to  influence  the  saddle  point  criterion,  we  obtain  a  formu¬ 
lation  that  includes  such  wave  phenomena  as  the  reflected  lateral  wave  field  in  the 
water  and  the  evanescent  field  in  the  bottom.  Each  saddle  point  of  a  given  plane 
wave  integral  corresponds  to  the  angle  of  an  eigenray  that  contributes  to  the  field, 
and  eigenray  characteristics  such  as  travel  time,  propagation  path,  and  reflection 
or  transmission  losses  are  derived.  The  simple  ray  interpretation  breaks  down  to 
a  certain  degree  when  the  saddle  point  represents  a  complex  eigenray  angle.  This 
occurs  when  the  field  is  composed  of  several  ray  contributions  merging  into  one, 
or  when  the  field  is  evanescent.  Although  saddle  point  approximations  to  the 
integral  for  the  field  are  usually  adequate,  a  technique  for  numerical  integration 
along  the  path  of  steepest  descent  must  be  used  in  certain  cases. 

The  prime  motivation  behind  developing  the  ray  method  for  the  single- 
interface  problem  is  to  apply  it  to  simple  shallow  water  ocean  environments.  In 
Chapter  3  we  extended  the  method  to  flat  and  sloping  waveguides  by  construct¬ 
ing  the  multiple  image  sources  that  account  for  the  multiple  ray  paths  between 
source  and  receiver  in  a  waveguide.  The  field  due  to  each  ray  path  is  constructed 
as  a  plane  wave  integral  in  a  manner  analogous  to  the  single-interface  problem. 
Each  integral  contains  multiple  reflection  coefficient  factors  (evaluated  at  differ¬ 
ent  angles  for  the  sloping  case),  corresponding  to  the  reflections  encountered  by 
the  ray  path  whose  contribution  the  integral  represents.  We  use  the  same  algo¬ 
rithm  to  evaluate  the  ray  fields  in  terms  of  saddle  points  ami  liteir  corresponding 
eigenrays,  and  we  obtain  the  total  field  in  the  waveguide  by  simply  summing  the 
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contributions  of  the  important  rays. 

We  have  evaluated  our  ray  model  algorithm  by  comparing  its  field  calcu¬ 
lations  with  those  produced  by  other  computer  models.  For  the  flat  oceanic  wave¬ 
guide,  the  ray  model  and  the  SAFARI  fast  field  program  predict  almost  identical 
fields.  Results  from  normal  mode  theory  are  also  compared,  and  we  find  that  if 
only  the  trapped  normal  modes  are  calculated  (as  is  usually  the  case),  the  mode 
results  can  be  poor  when  a  mode  is  near  cutoff.  We  find  good  agreement  between 
results  from  the  ray  model  and  the  SAFARI  model  even  in  the  very  low  frequency 
case  where  the  depth  of  the  water  is  only  half  of  the  acoustic  wavelength. 

For  the  sloping  waveguide,  or  wedge,  the  ray  model  has  been  applied  to 
benchmark  problems  proposed  in  December  1986  by  the  Underwater  Acoustics 
community  of  the  Acoustical  Society  of  America.  In  Chapter  4  we  demonstrated 
excellent  agreement  between  results  from  the  ray  model  and  a  model  based  on 
two-way  coupled  mode  theory.  The  independence  of  the  two  models’  approaches 
assures  us  that  their  common  solution  is  indeed  correct.  Since  there  are  always 
uncertainties  in  the  formulation  and  implementation  of  any  propagation  model, 
the  development  of  a  new  and  independent  technique  is  a  valuable  contribution  to 
the  modeling  effort  in  general.  For  example,  the  availability  of  the  coupled  mode 
results  was  much  appreciated  when  the  author  was  initially  implementing  the  ray 
algorithm  as  a  computer  program. 

Besides  constituting  an  independent  means  of  calculating  acoustic  fields 
in  shallow  water,  the  ray  model  provides  physical  insight  into  the  composition  of 
the  field  in  terms  of  eigenrays  and  their  characteristics.  In  Chapter  4  we  addressed 
the  issues  of  backscatter,  the  importance  of  lateral  waves,  and  the  importance  of 
partially  reflected  rays  by  analyzing  the  eigenrays  at  a  particular  receiver  position 
in  a  wedge. 

The  ray  algorithm  differs  from  most  numerical  models  in  that  it  can  find 
the  field  in  a  pointwise  manner,  as  opposed  to  “marching”  algorithms  that  must 
compute  the  field  at  all  points  between  a  given  source  and  receiver.  In  applications 
where  the  field  at  only  a  handful  of  points  is  of  interest,  the  pointwise  algorithm 
can  be  very  advantageous.  Another  unique  feature  of  the  ray  method  is  that  it 
is  more  efficient  at  high  frequencies  than  at  low  frequencies.  In  addition,  the  ray 
method  has  no  large  computer  memory  requirements  that  are  sensitive  to  the 
acoustic  frequency  being  modeled. 

Ray  modeling  of  propagation  in  a  three-dimensional  wedge  with  a  pen¬ 
etrable  bottom  was  described  in  Chapter  5.  In  deriving  the  plane  wave  integral 
for  each  ray  field  in  this  case,  we  find  that  a  double  integral  over  two  independent 
plane  wave  angles  must  be  evaluated.  As  a  consequence,  the  saddle  point  crite¬ 
rion  requires  partial  derivatives  with  respect  to  two  angular  variables  to  vanish. 
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A  more  complete  investigation  of  the  feasibility  of  solving  the  two-dimensional 
integral  remains  to  be  done,  but  an  approximate  method  that  would  be  straight¬ 
forward  to  implement  is  proposed.  Three-dimensional  propagation  problems  are 
of  great  interest  to  the  underwater  acoustics  community,  but  they  pose  serious 
problems  to  all  modeling  approaches.  Successful  application  of  the  ray  model  to 
even  the  simplest  three-dimensional  problem,  the  isovelocity  wedge,  would  provide 
a  valuable  tool  in  evaluating  other  modeling  methods. 

The  final  application  of  the  ray  method,  described  in  Chapter  6,  was  to 
model  a  directional  source  by  assigning  complex  spatial  coordinates  to  a  point 
source.  Although  a  point  source  in  complex  space  is  difficult  to  interpret  physi¬ 
cally,  the  mathematical  formulation  of  the  problem  results  in  branch  cut  surfaces 
in  real  space  that  can  be  viewed  as  containing  equivalent  sources.  The  reflected  and 
transmitted  fields  due  to  a  beam  incident  on  a  fluid-fluid  interface  are  computed 
using  the  same  techniques  as  for  the  point  source  problem.  The  only  modification 
required  is  that  certain  parameters  describing  the  source-receiver  geometry  must 
be  made  complex.  Several  examples  of  fields  due  to  beams  incident  near  the  crit¬ 
ical  angle  have  been  computed,  and  brief  analyses  have  been  given.  A  directional 
source  in  a  flat  or  sloping  waveguide  could  be  modeled  using  the  same  techniques. 

The  development  of  the  ray  methods  described  in  this  disseration  has 
focused  on  simple  environmental  models.  For  example,  the  ocean  water  has  been 
assumed  isovelocity,  interfaces  have  been  assumed  flat,  the  media  have  been  as¬ 
sumed  to  be  fluids,  and  ocean  bottoms  have  been  assumed  homogeneous.  As  a 
research  tool,  it  is  valuable  to  be  able  to  address  such  idealized  problems  with 
great  precision;  in  practice,  however,  it  is  desirable  to  describe  real  ocean  en¬ 
vironments  more  accurately  by  including  sound  velocity  gradients  in  the  media, 
interface  roughness,  and  shear  velocities  and  layering  in  the  ocean  bottom.  Fu¬ 
ture  research  may  involve  applying  some  of  the  concepts  and  numerical  methods 
developed  here  to  tackle  these  more  difficult  problems. 
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